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Abstract The brain's activity is characterized by the interaction of a very large num- 
ber of neurons that are strongly affected by noise. However, signals often arise at 
macroscopic scales integrating the effect of many neurons into a reliable pattern of 
activity. In order to study such large neuronal assemblies, one is often led to derive 
mean-field limits summarizing the effect of the interaction of a large number of neu- 
rons into an effective signal. Classical mean-field approaches consider the evolution 
of a deterministic variable, the mean activity, thus neglecting the stochastic nature 
of neural behavior. In this article, we build upon two recent approaches that include 
correlations and higher order moments in mean-field equations, and study how these 
stochastic effects influence the solutions of the mean-field equations, both in the limit 
of an infinite number of neurons and for large yet finite networks. We introduce a new 
model, the infinite model, which arises from both equations by a rescaling of the vari- 
ables and, which is invertible for finite-size networks, and hence, provides equivalent 
equations to those previously derived models. The study of this model allows us to 
understand qualitative behavior of such large-scale networks. We show that, though 
the solutions of the deterministic mean-field equation constitute uncorrelated solutions 
of the new mean-field equations, the stability properties of limit cycles are modified 
by the presence of correlations, and additional non-trivial behaviors including periodic 
orbits appear when there were none in the mean field. The origin of all these behav- 
iors is then explored in finite-size networks where interesting mesoscopic scale effects 
appear. This study leads us to show that the infinite-size system appears as a singular 



J. Touboul 

*NcuroMathComp Laboratory 

INRIA/ ENS Paris 

23 Avenue d'ltalic 

75013 Paris Tel.: +33-1 39 63 57 10 

Fax: +33-4 92 38 78 45 

E-mail: jonathan.touboul@sophia.inria.fr 

** G. Bard Ermentrout and Jonathan Touboul 

Department of Mathematics, 

University of Pittsburgh, 

Pittsburgh PA, 

USA. 



2 



limit of the network equations, and for any finite network, the system will differ from 
the infinite system. 

Keywords Neural Mass Equations; Dynamical Systems; Markov Process; Master Equa- 
tion; Moment equations; Bifurcations; Wilson and Cowan system. 

Introduction 

Cortical activity manifests highly complex behaviors which is often strongly character- 
ized by the presence of noise. Reliable responses to specific stimuli often arise at the 
level of population assemblies (cortical areas or cortical columns) featuring a very large 
number of neuronal cells presenting a highly nonlinear behavior and that are intercon- 
nected in a very intricate fashion. Understanding the global behavior of large-scale neu- 
ral assemblies has been a great endeavor in the past decades. Most models describing 
the emergent behavior arising from the interaction of neurons in large-scale networks 
have relied on continuum limits ever since the seminal work of Wilson and Cowan and 
Amari 2,3,53,54 . Such models tend to represent the activity of the network through 
a macroscopic variable, the population-averaged firing rate, that is generally assumed 
to be deterministic. Many analytical and numerical results and properties have been 
derived from these equations and related to cortical phenomena, for instance for the 
problem of spatio-temporal pattern formation in spatially extended models (see e.g. 

This approach tends to implicitly make the assumption that correlations in the 
activity do not modify the behavior of the system in large populations. However, the 
stochasticity of the system, and in particular the presence of correlations in the activity 
may considerably affect the dynamics in such nonlinear systems, and these can be 
even more important as the network size increases. These models therefore make the 
assumption that averaging effects counterbalance the prominent noisy aspect of in vivo 
firing, a feature that can dramatically affect finite-sized network activity, and, thus, 
that the emergent behavior of a cortical column is deterministic. 

However, increasingly many researchers now believe that the different intrinsic or 
extrinsic noise sources are part of the neuronal signal, and rather than a pure disturb- 
ing effect related to the intrinsically noisy biological substrate of the neural system, 
they suggest that noise conveys information that can be an important principle of 
brain function 46 . At the level of a single cell, various studies have shown that the 
firing statistics are highly stochastic with probability distributions close to Poisson 
distributions [47], and several computational studies confirmed the stochastic nature 
of single-cells firings 15,43,50,51 . How variability at the single neuron level affects 
dynamics of cortical networks is less established so far. Theoretically, the interaction 
of a large number of neurons that fire stochastic spike trains can naturally produce 
correlations in the firing activity of the population considered. For instance power-laws 
in the scaling of avalanche-size distributions have been studied both via models and 
experiments 7,8, 37, 49j . In these regimes the randomness plays a central role. 

A different approach has been to study regimes where the activity is uncorrelated. 
A number of computational studies on the integrate-and-fire neuron showed that under 
certain conditions neurons in large assemblies end up firing asynchronously, produc- 
ing null correlations 1,4,14 . In these regimes, the correlations in the firing activity 
decrease towards zero in the limit where the number of neurons tends to infinity. The 
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emergent global activity of the population in this limit is deterministic, and evolves 
according to a mean-field firing rate equation. However, these states only exist in the 
limit where the number of neurons is infinite, and such asynchronous states do not 
necessarily exist in cortical areas or in computational models, and even when they 
exist, they are not necessarily stable. This raises the question of how the finiteness of 
the number of neurons impacts the behavior of asynchronous states, and how to study 
non-asynchronous behaviors. The study of finite-size effects in the case of asynchronous 
states is generally not reduced to the study of mean firing rates and can include higher 
order moments of firing activity ID 22 b! . In order to go beyond asynchronous states 
and take into account the stochastic nature of the firing and how this activity scales 
as the network size increases, different approaches have been developed, such as the 
population density method and related approaches [18]. Most of these approaches in- 
volve expansions in terms of the moments of the resulting random variable, and the 
moment hierarchy needs to be truncated which is a quite hard task that can raise a 
number of technical issues (see e.g. [38]). A recent approach overcomes this truncation 
difficulty by considering the mean-field behavior as a stochastic process and therefore 
treating the rigorous limit equation of the interaction of many stochastic neurons as a 
fixed point equation in the space of stochastic processes, making use of the powerful 
tools of stochastic limit theorems and large deviation techniques (see [27] ) . 

In order to study the effect of the stochastic nature of the firing in large network, 
many authors strived to introduce randomness in a tractable form. The models we 
will be interested in are based on the definition of a Markov chain governing the fir- 
ing dynamics of the neurons in the network, where the transition probability satisfies 
a differential equation called the master equation. Seminal works of the application 
of such modeling for neuroscience date back to the early 90s with the work of Ohira 
and Cowan [42], and have been progressively developed by different authors and today 
constitute a model of choice 16 ,17,10, 22] . The present manuscript is based on the anal- 
ysis of two recently developed stochastic models of the customary Wilson and Cowan 
equations. The global resulting activity is then readily derived from the activity of all 
neurons in the network. The two models of interest here were instantiated specifically 
in order to be compared to standard Wilson and Cowan rate equations, and provide 
tractable stochastic extensions of the purely deterministic standard neural-field models 
|17II10] . The two approaches differ in the choice of the transition probability and on the 
variables considered, choices mainly based on the fact that they address two different 
regimes of network activity: Buice, Cowan and Chow address the highly correlated 
Poisson-like activity and Bressloff, the asynchronous states. The two approaches also 
differ in the way the equations are treated: Buice, Cowan and collaborators in [16, 17 
derive moment equations by carrying out a loop expansion of a path integral represen- 
tation of their master equation whose truncation is justified under the assumption that 
the network operates in a Poisson-like regime. Bressloff rigorously justifies his moment 
truncation by exhibiting an expansion in powers of a small parameter 1/N where N 
corresponds to the size of the network. These different points of view lead the authors 
to define different variables and they obtain closely related yet different equations on 
the variables they define, derive equations for moments from the master equation which 
they then simplify and truncate to finite order. 

The present paper is organized as follows: in the first section we describe the models 
we use. Both simplified models introduced involve coupling between mean firing rate 
and correlations between firing times, with finite-size correlations. In the second section, 
we are interested in the non-zero correlation solutions of the mean-field equations, and, 
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in particular, focus on the qualitative similarities and differences between the solutions 
of these stochastic equations compared to the customary Wilson and Cowan behaviors. 
In the third section, we turn to study the finite-size effects, namely the qualitative 
differences between the solutions of mean-field limits and Wilson and Cowan system on 
one hand, and the behaviors of large but finite networks on the other hand. These finite- 
size solutions are compared back to the initial Markovian model and the differential 
equations they produced. 

1 Markovian firing models 

In this section, we introduce and describe the master equation formalism used by 
Buice, Cowan and collaborators [16|ll7j (hereafter BCC) and Bressloff [10] and extend 
these approaches to take into account different populations. Both models have been 
previously derived, so we quickly present their derivation and the resulting equations in 
this section, as our main purpose is to derive from these equations our model of interest, 
the infinite-size model, and to mathematically analyze it. Details on the derivation of 
these equations and on their theoretical justification can be found in [16II17IITU] . 

Throughout manuscript, we will be interested in a network structured into M 
homogeneous neural populations. Each population i £ {1, . . . , M} is composed of TVj 
neurons, and the total number of neurons is denoted by N (i.e. N = J2t=i N i)- We 
are interested in the behavior of the network when the number of neurons TV tends 
to infinity. In this limit, we assume that the proportion of neurons belonging to each 
population are non-trivial, i.e. lirn/v->.oo N{/N = A, £ (0, 1). Each neuron of population 
i interacts with all the neurons of the network, with an efficiency denoted Wij . These 
weights are assumed to scale as l/Nj, so that a given activity at the level of population 
i will produce a bounded activity at the level of population j. We therefore define the 
effective interconnection weights Wij such that Wij = Wij/Nj. Each population is 
assumed to receive an input firing rate h(t) considered deterministic and constant in 
the rest of the paper. 

1.1 Deterministic Wilson and Cowan equations 

The Wilson and Cowan (WC) model is based on the assumption that the synaptic input 
current to each population is a function of the firing-rate of the pre-synaptic population, 
and that the contributions of the different populations are linearly summed to produce 
the post-synaptic population firing rate. It assumes moreover that the population- 
averaged firing rates Vi(t) are deterministic function of time that are governed by the 
equations: 



where fi is a function transforming an incoming current into an output firing rate 
and typically has a sigmoidal shape, ccj are the relaxation rates corresponding to the 
natural inactivation of each population when they receive no input. 

This purely deterministic approach, often used for spatially extended neural net- 
works, has proved efficient to model a large number of cortical phenomena, but fails to 
reproduce the stochastic behaviors that can appear at the level of population activity. 
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We now describe the stochastic framework used by ll6H71iT0] to extend this framework 
in order to take into account the stochastic nature of the firing. 



1.2 Markovian Framework 

In the Markovian approach yielding BCC and Bressloff models, each neuron can be 
either quiescent or active (active meaning in the process of firing an action potential), 
and the state of the network is then described by the variable n £ N M where rii{i) is 
the total number of active neurons in population i. This chain is modeled as a one- 
step discrete-time Markov process, i.e. the only possible transitions of this chain are 
assumed to be rij — > m ± 1 for i G {f , . . . , M}. Each active neuron of population 
i returns to a quiescent state with constant probabilistic rate on and each quiescent 
neuron become active with a state-dependent rate Fi(n) depending on the inputs the 
neuron receives. The probability of the network to be in a state n at time t therefore 
satisfies a differential equation called the master equation: 

M 

i=l 

+ Fi (n^)P(n^ ,t)-Fi (n) P(n, *)] (2) 

From this equation the authors derive equations on the moments of the Markov chain. 
We propose in appendix|A]an alternative derivation of the moment equations based on 
a Langevin approximation of the rescaled Markov chain and on the use of Rodriguez 
and Tuckwell expansions [44 . 

In the present paper we are particularly interested in a variable that would corre- 
spond in the Markov setting to the proportion of active neurons in each population: 
Pi(t) — n,-(t)/iVj for i = 1, . . . , N, which is a priori a well-behaved bounded stochastic 
process that takes values in the interval [0, 1]. Having defined this variable, we will be 
in a position to rescale BCC and Bressloff variables, and obtain from their derivation 
the equation governing this variable. In doing so, as further illustrated below, we will 
not need to derive new moment equations from the master equation, but it will ap- 
pear as a change of variable in the set of ordinary differential equations these authors 
previously obtained. In this view, the rescaled mean firing rate is defined as the instan- 
taneous averaged proportion of active neurons in each population fj(t) = (nj(t))/iVj, 
and rescaled correlations as the average value of the correlations of the proportion of 
active neurons in each population Cjj = {n,i(t)nj(i)/NiNj), the two first moments of 
the process p(t). 

The equations for the moments involve averaged values of (Fi(n)) which, since the 
activation is not polynomial, involve all the moments of the variable n, and therefore 
need to be truncated in order to obtain a closed set of equations. The choice of the 
function Fi is so far unspecified and will be defined in order to recover the WC equations 
in the mean-field uncorrelated limit. 

With the aim of specifically investigating (uncorrelated) asynchronous states, Bressloff 
truncates the moment series expansion up to second order, and considers the evolution 
of two coupled variables: the mean firing-rate Vi and riij = NjCij a correlation vari- 
able differing from our variable in that it is scaled by Nj . After Taylor expansion and 
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a Van-Kampen system-size expansion, a specific choice of the activation function Fi is 
made so that one ends up with the following set of ordinary differential equations: 

^ = -Oi m + fi(si) + fi( s i) Hk.l w ik w il n kl 

' ~3t~ =[ a iVi + fi(si)]Sij - (ai + aj)mj (3) 
+ J2klf'i( S i) W ik n k 3 + fj(sj)w jk n ki ] 

where Vi corresponds to the firing rate of population i: Viif) = {rii(t))/N and nij 
are defined as the correlations (rii(t)nj(t)}/N. In the original formula, the number of 
neurons of each population is assumed to be equal to TV. This quantity is bounded by 
application of Kurtz theorem [34], and therefore the term i fl'( s i) X^fc I w ik w il n kl 
tends to zero as TV tends to infinity, implying that the equation of the mean becomes 
independent of the correlations riy. Since we are interested in the statistics of the 
proportions of active neurons in each population, we are in a position to derive from 
the original Bressloff equations (J3J) the system of equations governing our variables Vi 
and Cij formally, and obtain in a multi-population case with distinct population sizes: 

= -CiM + fi(si) + \ fi(Si) Y^ k ,l w lkWilCkl 
' ^dT~ = Wi\ a i V i + fi( s i)] S ij ~ ( a i + a ]) C ij ( 4 ) 
+ Efe [f'i(. s i) w ik Cfej + fj (Sj)wjk Chi] 

In both equations |3| and Q, we denoted by the total instantaneous current 

received by population i at time t: 

M 

s i {t) = y £ j w ij v j {t) + I i {t). (5) 

We note that the system Q, though derived from Bressloff 's master equation and 
formalism, is equivalent to the original Bressloff model |3| as long as the total number 
of neurons N is finite, since it was derived through an invertible change of variable 
which is invertible as long as TV is finite, but not invertible in the limit TV — > oo. 

Buice, Cowan and collaborators [16,17 (BCC case), interested in studying the 
Poisson-like firing modes of the network, transform the moment equations to derive 
equations of normal ordered cumulants measuring the deviations of the moments of 
the variable n from pure Poisson statistics. The first normal ordered cumulant is equal 
to (nj), the second ordered cumulant to Cjj = — Mj/iVj5y. The initial approach 
of Buice, Cowan and collaborators is not a finite-size expansion per se in the general 
case, hence the parameter according to which the expansion is performed is no longer 
1/TV as in the case of Bressloff 's Van Kampen expansion, but consists of a multiple 
time scale expansion where the small parameter is the decay time of normal order 
cumulants (which applies far from bifurcation points). However, in the fully connected 
case we are interested in in the present manuscript with weights scaling as 1/TV, the 
expansion provided Buice, Cowan and collaborators method can be reduced to a finite- 
size expansion, and can be compared to Bressloff's case. In terms of these variables 
after some approximations and moment truncation to the second-ordered cumulant, 
they end up with the following set of ordinary differential equations (this step involves 
an instantiation of the activation functions Fj different from the Bressloff case): 

= -am + fi(si) + \ fi'{Si) Yjj k w ij w ikCjk 

d 

< = -(a, t + aj)cij + f[{s,i) J2 k w lk c k j + fj{sj) J2 k w jkCki+ (6) 

+ Wjfi( s i) w ij v i + JV7-/j( S j) W H V i 



7 



Remark 1 We note that the rescaling does not affect the form of the equations obtained 
in [T7] . These equations are valid under our assumptions that the synaptic weights WJy 
scale as 1/N in a fully connected network. We emphasize the fact that both expansions 
can be rigorously derived far away from bifurcation points. In the manuscript, we will 
however be interested in the bifurcations of the dynamical systems given by Q and 
Q. Indeed, these bifurcations produce changes in the number and in the stability of 
attractors (fixed point and limit cycle) which will be visible in non-trivial parameter 
intervals, and therefore will affect the behaviors of the system far away from bifurcations 
and will give us indications of the number and stability of attractors in those parameter 
regions. 

These models formally appear to very similar. However, it is important to note 
that the these two approaches are different: the master equation corresponds to differ- 
ent transition rates (different choices of the functions Fi as stated), and the variables 
considered differ. The two equations Q and (JsJ) constitute our starting point for a 
mathematical exploration of the solution, and in the rest of this paper, we mathe- 
matically analyze these equations. The reader nevertheless needs to keep in mind the 
differences between the two approaches for interpreting the results. 



1.3 The infinite-size model 

The equations Q and |6]| are a priori different. First of all, they do not deal with 
the same quantities: equation Q couples the mean firing-rate and the correlations 
while equation (JsJ) the mean firing-rate and the first order cumulant. However, in the 
limit N — > oo, the first-order cumulants, Cjj = Cij — Vi/Ni5u, are simply equal to 
the correlation Cy. Moreover, in the rescaled models we introduced, we observe that 
when the number of neurons tends to infinity, both Bressloff Q and BCC Q models 
converge to the same equations. These equations will be referred to as the infinite size 
model and are given by the equations: 

^ = -OLiUi + fcfa) + lfi(si) J2 k i w ik^ l iA kl 
' ^ + _ (7) 

+ J2klfi( S i) w ikAkj + f'j(Sj)w jk A ki ] 

We recall that although the interpretation of the variable A differs in the general case, 
it is the not the case of the infinite model: in BCC model Aij — are the second 
ordered cumulant, while in Bressloff model A^j = Cij is the correlation in the firing 
activity, and in the infinite-size limit, the normal ordered cumulant is equal to the 
correlation. However, when N becomes finite, the reader needs to bear in mind that 
in the finite-size unfolding of the infinite-size equations in BCC case, A represents 
the deviation of spike statistics from a Poisson Process and in Bressloff model the 
covariance of the firing activity. When letting N be finite, the two models unfold the 
behavior of the system in a different fashion. We also note that if the variable A(t) is 
equal to zero, then the mean-firing rates Vi(t) satisfy the WC equations. In that view, 
Bressloff and BCC models are generalizations of WC system that take into account the 
correlations in the firing. 

In this paper, we will first be interested in studying the mean-field limit solutions of 
the infinite size equations |7|), and second in unfolding these mean-field behaviors when 
the number of neurons is large but finite. Our study will particularly focus on two main 
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aspects: (i) how the stochastic nature of the firings affect the observed behaviors at the 
macroscopic level through the interplay of firing activity and the correlations, and (ii) 
how the finiteness of the number of neurons in the network disturbs these behaviors. 

We also note that in view of Kurtz' theorem [34], our variable in the infinite-size 
system vanishes. We formally consider these equations independently of this restriction. 
In other words, the infinite-size system is obviously not equivalent to the WC system, 
but BCC and Bressloff systems should be equivalent, in this limit, to WC system. The 
study of the infinite-size system will allow us to find what we call correlation-induced be- 
haviors, which are qualitative distinctions between the behaviors of the solutions of the 
infinite-size system and of Wilson and Cowan system. If these behaviors are unfolded 
non-trivially into solutions of the finite-size equations, these will evidence additional 
behaviors of the system presented by the finite-size Bressloff and BCC equations that 
do not correspond to solutions of Wilson and Cowan system. We will address these 
questions in in section [3] 

Besides our theoretical analysis on the equations, we will also go back to the initial 
Markovian model and compare (with Monte Carlo simulations) the stochastic behaviors 
it presents to the solutions of WC, BCC and Bressloff models. 



2 Influence of the correlations in the infinite-size system 

We have shown that in the limit where number of neurons is infinite, both Bressloff 
and BCC models, in our particular rescaling, converge to the infinite-size model given 
by equations |7| where the mean-firing rate is coupled to the correlations in the firing 
activity. When the correlations are equal to zero, the equation for the mean-firing rate 
is uncoupled from the one for the correlations A and is identical to the WC system. The 
question we address in this section is how the interaction between the mean activity and 
the correlations modifies the behavior of the global activity compared to WC system 
given by equations |TJ. We address the following two questions: (i) are the solutions of 
Wilson and Cowan equations also solutions of the full systems including correlatons, 
and if the answer is yes, are the stability properties of these solutions conserved? and 
(ii) Do the correlations yield other behaviors? 



2.1 Wilson and Cowan's Solution Solve the Infinite-Size System 

In equation Q, it is easy to see that zero correlation: A(t) = is always a fixed 
point of the correlation equations no matter what mean firing rates v(t) are. For zero 
correlations, the equations for the mean firing rate v are reduced to the classical WC 
equations. Therefore, WC solutions define solutions for the infinite-size system with 
zero correlations. Let us now study the stability of the solutions defined by these 
uncorrelated WC behaviors in the infinite size system. 

The infinite-size system involves a standard M-dimensional differential equation 
on the mean firing rates v(t) coupled to a matrix differential equation. In order to 
use customary approaches for dynamical systems, we transform the correlation matrix 
into a M 2 dimensional vector and express the equations on A in this new format using 
Kronecker products from linear algebra. To this end, we start by defining the function 
Vect transforming a M x M matrix into a M 2 -dimensional column vector, as defined 
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in 1411 



Vect 



MxM 



^ [Xn,. 



,X M1 ,Xi2(t), . . . ,X M2 (t), . . . X 1M (t), . . . ,X MM (t)] 

mxn and B e n rxs 



Let us now denote by <g> the Kronecker product defined for A 6 1R 
as the (m r) x (n s) matrix: 



>B 



( a\\B ayiB 

a 2\B CI22B 



a ln B \ 
a 2n B 



\ a m \B a TO 2-B ■ • ■ a mn B J 



Some definitions and identities in the field of Kronecker products are reviewed in ap- 
pendix [B] and a more complete discussion can be found in and references therein. 
We recall here for the sake of completeness some well-known relationship that will be 
useful. Let A, B, D,G,X g ]R MxM 5 let I M be the MxM identity matrix and A ■ B 
or AB denote the standard matrix product. We have: 

' Vect(AXB) = (B T ® A)Vect(X) 

Ae A = A® I M + I M ® A . (8) 

(A <g> B) ■ (D <g> G) = (A ■ D) x (B ■ G) 

The relationship © is called Kronecker sum. In this framework, we have the following 
identity: 

Lemma 1 Let V A (t) = Yect(A(t)) be the column vectorization of the matrix A(t), and 
A(x) the matrix defined by A(x) — —a + F(x) where a is the diagonal matrix with 
element am — on and F the matrix of general element {F)ij = (/iQ^feLi w ik x k) w ij)- 
The variable V A (t) satisfies the differential equation in IR^ 
dV A {i) 



dt 



= (A(v(t))(BA(v(t)))VA(t). 



(9) 



Proof The differential equation governing the evolution of the coefficients of the corre- 
lation matrix A(t) of the infinite size system |7| can be easily reordered into: 

M 



dAi. 
dt 



(A ik (u(t))A kj (t) + A jk (v(t))Ai 



fe=i 



which, through straightforward linear algebra manipulations, can be written as: 

^ = A(v{t)) ■ A(t) + A(t) ■ A{v{t)) T . 

The linear operator X H> A(u(t)) ■ X + X ■ A(u(t)) T for X £ ]R MxM can be written 
in terms of Kronecker product of matrix on the vectorized version V A (t) of the matrix 
A{t) using the relationship given in (JsJ and we have: 

Vect(A(v{t)) ■ A(t) + A(t) ■ A(v(t)) T ) = Vect (,40(*)) ■ A(t) ■ I M ) +Vect(l M ■ A(t) ■ A(u(t)f) 

= {i^ ® A(u(t)) + A(v(t)) ® I M ) ■ V A (t) 
= (A{u(t)) © A{u(t)))-V A (t). 
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Now that we have defined a convenient framework to study the infinite size equa- 
tions, we are in position to study the stability of fixed points and limit cycles of the 
WC system as solutions with zero correlations of the infinite size system. We start by 
addressing the problem in the one population model (M = 1), which allows straight- 
forward calculations and will be helpful in our specific study of the one population 
infinite-size equation. 

Proposition 1 In the one population case, any fixed point of Wilson and Cowan equa- 
tion is a zero- correlation fixed point of the infinite-size system, and with the same 
stability properties. 

Proof In the one population case, equations (JsJ) constitute a planar dynamical system 
that reads: 



where s = w v + 1. 

As stated, the solution A — is a fixed point of the correlation equation, and in 
that case, the activity satisfies the one-population WC equation: 



Let v* be a fixed point of this system. The Jacobian matrix of the infinite-size system 
at the fixed point (u* , A = 0) reads: 



As we can obviously see, the eigenvalues of this equilibrium are Ai := (— ct + w f (wv* + 
I)) and A2 = 2Ai, and therefore the stability of this fixed point is the same as the 
stability of the related fixed point in WC's system. 

We therefore conclude that any stable solution of the one dimensional WC equa- 
tion provides a zero correlation (A = 0) stable fixed point of the infinite-size system, 
and we observe that no stable solution of WC system is destabilized by the presence 
of correlations, and no unstable solution of the Wilson and Cowan equation will be 
stabilized. Therefore, zero-correlation fixed points of the infinite size system exists if 
and only if Wilson and Cowan equation has a fixed point, and the related fixed point 
of the infinite-size system has the same stability as WC's fixed point. 

We now turn to demonstrating the related properties in arbitrary dimensions M, 
for fixed points and cycles. 

Theorem 1 Solutions of the WC system provide solutions of the infinite-size system 
with zero correlations A — 0. The stability of fixed points and limit cycles of the infinite 
size system depend on the stability of the solution for the WC system as follows: 

i) A fixed point of the infinite-size equation with null correlations is stable if and 
only if the value of related mean-firing rate is a stable fixed point of WC system; 




(10) 



IE 



av + f(w v + I) 




li) A cycle of the infinite size system with zero correlations (v(t),A = 0) is expo- 
nentially unstable if and only if v(t) is an unstable cycle of WC system. Stable cycles 
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v(t) of WC's system define cycles of the infinite-size system with zero correlations with 
neutral linear stability. The infinite size system does not present any stable cycle with 
null correlations. 

Proof In order to prove the theorem, it is convenient to introduce the vector fields 
F v : IR M H> R M , Gu{v) : 1R M i — ^ i? MxM " and F A ■ H M h-> jrA^xA/ 2 that gQvern thg 
dynamics of the infinite size system: 

% = F v {y) + Gv{v)V A {t) 
= F A {v)V A (t) 

We have F v (y) = -a + [f i {s l ),i = 1 . . . M] T , F A (v) = A(v) © A(v) where A(-) is 
defined in lemma [l] and Gy{v) is the M 2 x M 2 matrix such that: 



Gu(v)V A (t) = [\ fi(8i) w ik w u A kl ,i = 1 . . . Mf 

k,l 



Let {y{t), A = 0) be a solution of the infinite-size system. The Jacobian matrix of 
the system evaluated on this solution has the form: 



J(u(t),0) 



d v {Fu) + d v {G v )V A \ G v {u)\ 

\u(t),V A =0 \v(t),V A =0 
du{F A )V A \ d VA (F A )\ 



= ( du(F v )(u(t)) G v {u{t)) \ 
V d VA (F A )(u(t)) J 

where dx f for a multidimensional vector X and a multidimensional function / denotes 
the differential of / with respect to X. Since F v is exactly the vector field associated 
with WC system, d v F v is an M x M block corresponding to Jacobian matrix of WC 
system at the point v, which is equal to the matrix A(y) introduced in lemmajl] 

Now that the Jacobian matrix is identified, we prove the assertions of the theorem: 

i) Let us consider a fixed point {v, A = 0) of the infinite size system. Necessarily, v 
is a fixed point of Wilson and Cowan equation. The stability of this fixed point depends 
on the spectrum of the Jacobian matrix J of the system at this point. Since the Jacobian 
matrix of the infinite-size system is block diagonal, its spectrum is therefore composed 
of the eigenvalues each diagonal block d v (Fv) — A(v) and dy A (F A )(v(t)). The first 
block, A{v), has exactly the eigenvalues of Wilson and Cowan system at its fixed point 
v. These eigenvalues are denoted {Xi, i = 1 . . . M}. The second diagonal block of the 
Jacobian matrix dy A (F A ) is equal to the Kronecker sum A(v) © A(y) by application of 
lemma[l] The eigenvalues of a Kronecker sum are known to be all possible pairwise sums 
of all the eigenvalues of A(y), viz. (Aj + Xj; (i,j) G {1, . . . , M} 2 ) (see an elementary 
proof of this fact in proposition [2] of appendix [B|, hence the spectrum of the Jacobian 
matrix of the infinite size system is composed of the eigenvalues: 

{\ i ,i = l...M}U{\ i + \ j ,{i,j) 6 {1,...,M} 2 }. 

These eigenvalues depend in a simple fashion on the eigenvalues of the WC Jacobian 
matrix at the fixed point v, and it is easy to show that the fixed point (y, 0) in the the 
infinite system is: 



12 



— exponentially stable if and only if the v is an exponentially stable fixed point of 
WC system, i.e. if and only if all the eigenvalues \ have a strictly negative real 
part. 

— exponentially unstable if and only if the v is an exponentially unstable fixed point 
of WC system, i.e. if and only if there exists an eigenvalue with strictly positive 
real part. 

— neutrally stable if and only if v is neutrally stable for WC system, i.e. if and only 
if all eigenvalues have non-positive real part and at least one eigenvalue have a null 
real part. 

In summary, the stability of the zero-correlation fixed point (v, A = 0) is exactly the 
same as the stability of v as a fixed point of WC system. 

ii) Let us now address the case of cycles will null correlations. To this end, let us 
consider that (v(t),A(t) = 0) is a periodic orbit of the infinite-size system (i.e. v{t) is 
a periodic orbit of a WC system). Let us denote by T > the period of this cycle. We 
show that the null correlations A = is an exponentially unstable solution if the cycle 
v(t) is exponentially unstable as a solution of WC system, and neutrally stable if the 
cycle v(t) is exponentially stable or neutrally stable as a solution of Wilson and Cowan 
system, which proves the theorem. 

Given that the cycle v(t) is known, the correlations A(t) satisfy a linear equation: 



A{v(t))®A{v{t)))V A {t) 



dV A (t) 
dt 

and the time-dependent matrix \ A{y{t)) © A(u(t))J is T-periodic. A basic result of 
Floquet theory implies that the stability of the solution V A = depends on the eigen- 
values of the resolvant of the system at time T, namely the matrix <?(T) in K x 
where $■(•) is defined by 

=(A(v(t))®A( V (t)))#(t) 



HO) = hn 

More precisely, the null fixed point is exponentially stable if and only if all eigenvalues 
of if'(T) ( the multipliers) are of modulus strictly smaller than 1, neutrally stable if 
all have modulus smaller or than equal to 1 with at least one multiplier with modulus 
equal to 1, and exponentially unstable if there exists a multiplier with modulus larger 
than 1. We therefore need to characterize the eigenvalues of ^(T) in order to conclude 
on the stability of the solution in the infinite-size system. To this end, let us introduce 
$ (t) the resolvent of WC system, i.e. the unique solution of the fundamental equation: 

<P(0) = I M 

We prove in theorem [3] of appendix [B| that &(£) = <P(t) ® <P(t). Let us denote by 
{lH, i = 1, ...,M} the eigenvalues of <5(T). Because ty(T) is the Kronecker square 
of <P(T), its eigenvalues are are all pairwise products of the eigenvalues of $(T), i.e. 
fij, i,j = 1...M}, as shown in proposition [2] of appendix [B| 

If the cycle v(t) is exponentially unstable as a solution of WC system, then Floquet 
theory implies that there necessarily exists a multiplier of <P(T), m, such that \m\ > 1. 
Therefore the resolvant ^(T) has a multiplier equal to $ which has modulus equal 
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to \m\ 2 > 1 and therefore this cycle is exponentially unstable as a solution of the 
infinite-size system. 

If u(t) is exponentially stable or neutrally stable as a solution of WC system, then a 
basic result of Floquet theory states that the resolvant if the linearized system evaluated 
on the cycle at time T, 4>(T), has at least a multiplier equal to 1 and all other multipliers 
with multipliers of modulus smaller than 1. We now provide an elementary proof of 
the existence of the 1 multiplier: let £(t) = F v {v(t)). Then since v(t) is T-periodic, so 
is and furthermore we have: 

dg«) = dfy^)) = {dvFv){v{t)) d^o = {dvFv)[v{t))Fv{v{t)) = A (v(t))m 

Moreover, since v(t) is not a fixed point, F(i/(Q)) 7^ 0. We have shown £(t) is solution of 
the linearized equation and hence using its T-periodicity, we have £(T) = #(T)£(0) = 
5(0), which proves that £(0) is eigenvector of $(T) associated with the eigenvalue 1. 
Therefore, the resolvent tf'(T) necessarily has an eigenvalue equal to one (associated 
with the eigenvector Vect(£(0)£(0) T )), which implies that the null fixed point has a 
neutral linear stability as a solution of the infinite-size system. We directly conclude 
from this result that the infinite size system does not features any exponentially stable 
cycle, which ends the proof. 

We conclude that the solutions of WC system always provide solutions of the 
infinite-size system, and that all fixed points with null correlations in the infinite-size 
system are fixed points of WC system, with the same stability properties. However, 
cycles of WC system lose exponential stability in the infinite-size system. Note that 
this does not necessarily implies that these cycles become unstable, and a nonlinear 
stability analysis is necessary. But it implies that the transient phase of convergence 
towards the cycle or of repulsion from the cycle will be not be exponential. 



2.2 Correlation-induced behaviors 

In the previous section, we just proved that all the fixed points of the WC system 
are solutions of uncorrelated activity in the infinite size system. We now investigate 
the existence of new solutions induced by the presence of non-null correlations in the 
infinite system, that are therefore not solutions of WC system. Such solutions will be 
referred to as correlation-induced behaviors. The general resolution of this problem is 
quite hard. For this reason, we provide a rigorous analysis of the one population case, 
and treat two-population case through numerical analysis and simulations. 

2.2.1 One population case 

We have seen in the one population case that any fixed point of WC system was solution 
of the infinite size system with zero correlations. Moreover, there is no possible cycle in 
the Wilson and Cowan equations in one dimension. We investigate now the existence 
of cycles or additional fixed points in the infinite size system with one population. 
We recall that acceptable solutions in the one population case necessarily have non- 
negative firing rate v and A. Indeed, in the case of the infinite-size model arising in 
the limit of Bressloff 's rescaled model, the correlation A needs to be positive as a limit 
of positive quantities (and as the covariance of the Markov process), and in BCC case 
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A is equal to the limit of the first order cumulant which differs from the correlation up 
to the coefficient ~v/N that vanishes in the infinite-size limit. 

Theorem 2 In the one population infinite size system, there does not exist any ac- 
ceptable stable fixed point with strictly positive correlation A. 



Proof Let us assume that A^Q. The fixed point equation of the infinite-size system 
-a + / (w v + I)w — 



-ai> + f(wu + I) + \f"{wv + I)w 2 A = 



2 . .. (") 



The first equation of the system (111 is independent of A and fixes possible values of 
v. Assuming that this equation has a solution v* and denoting s* = w v* + I, the 
Jacobian matrix of the system at v — v* reads, using the property that v* solves the 



first equation of (111 



1 rill I *\ 3 A 1 f" I *\ 2\ 

2f"(s*)w 2 A J { ' 

Such values of v* necessarily exist for some parameters, and since the function / is 
assumed to be sigmoidal, it can have multiple solutions. Moreover, because of the 
particular shape of the sigmoidal function /, the second derivative of / can vanish at 
a point corresponding to the inflection point of the sigmoid. 

In a very particular case (for precisely tuned values of the parameters), one can 
hence have / (w v* +1) — 0. In that case, fixed points only exist if — a+f (wu* + I)w = 
0, if this condition is satisfied, any value of the correlation A provides a fixed point of the 
system. This fixed point is never exponentially stable, since in that case the Jacobian 
matrix has a zero eigenvalue (obvious when using the expression of the Jacobian matrix 



1 12 1 and using the fact that f"(s*) = 0). 

In the general case where f"(s*) 7^ 0, the system has a fixed point with non-null 
correlations 

w 2 f"(s* 

We therefore need to check whether if this fixed point is stable and acceptable (i.e. 
mm(v* , A*) > 0). The Jacobian matrix at this point has the expression given by 
(12 1. Its determinant is equal to —(/"(«*)) w A A* and has an opposite sign to A* 



Therefore, for acceptable solutions with A > 0, the determinant of the Jacobian 
matrix is strictly negative. We conclude that any fixed point of the infinite-size system 
with A > are saddle fixed points. 



We conclude from this theorem that the system does not features any acceptable 
stable fixed point with non-zero correlations. We can now easily conclude that there are 
no acceptable cycles. Consider, first, a limit cycle such that A is strictly positive. Such 
a limit cycle must surround one or more fixed points whose Poincare index must sum 
to +1 5 . However, every fixed point in the postive orthant is a saddle point and saddle 
points always have indices of -1. Thus, no limit cycle can have a strictly positive A. 
Any limit cycle must therefore be tangent to the A = line and the point of tangency 
cannot contain a fixed point, so, again, the summed index of points in the limit cycle 
cannot be +1. 
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In conclusion, there are no acceptable limit cycles and all acceptable fixed points 
are saddles. Thus, the only stable behavior are the uncorrelated stationary states of 
the scalar WC equation. This situation will be quite different in the multi-population 
case, as we now develop. 

2.2.2 Correlation-induced behaviors in multi-populations networks 

We now turn to study multi-population networks with the aim of identifying solutions 
of the infinite-size system that qualitatively differ from WC equation. The situation in 
the multidimensional case will be quite different, since we have proven in theorem [T] 
that cycles of WC system lost exponential stability in the infinite-size system, allowing 
the appearance of distinct transient and/or asymptotic behaviors. 

To identify such correlation-induced behavior we numerically study two particular 
two population networks. The first model (Model I) is a two-population network in- 
cluding an excitatory and an inhibitory population known to produce oscillations, and 
the second example (Model II) builds on a famous model of binocular rivalry presenting 
bistability of fixed points. 

Perturbation of a WC cycle We proved in theorem [I] that any cycle of WC model lost 
exponential stability in the infinite-size system. This loss of exponential stability does 
not necessarily mean that the cycle loses stability, and this property depends on the 
higher order terms of the nonlinear equation, but even if the cycle keeps stability, the 
transient convergence phase towards the cycle will be dramatically slowed. 

In this section, we choose to study a two-populations network featuring an ex- 
citatory population interconnected with an inhibitory one. Because of the symmetry 
Aij = Aij, the system is of dimension 5 (the two mean firing-rates and three correla- 
tion variables). To fix ideas, we denote by 1 (resp. 2) the excitatory (resp. inhibitory) 
populations. We choose the same activation function for both populations equal to 
f(x) = 1/(1 + exp(— a;)) and the same inactivation rate a = 1, and are interconnected 
through the following connectivity matrix: 



Each population receives different input currents I\ = —0.5 and I2 = —5. These param- 
eters define a model called Model I in the sequel. For these functions and parameters, 
the WC system features a cycle. When taking the input to the excitatory population 
I\ small, the system features a single stable fixed point, that loses stability as the input 
I\ increases through a supercritical Hopf bifurcation generating a family of stable limit 
cycles, that disappear through a homoclinic bifurcation (see Figure [TJ a)), branching 
onto a high-state stable fixed point. Let us now analyze how this picture is modified 
by the presence of correlations. 

First of all, from theorem [l] it is clear that the parameter point where WC sys- 
tem undergoes the supercritical Hopf bifurcation will be a very degenerate bifurcation 
point for the infinite-size system. Indeed, the Hopf bifurcation is characterized by the 
presence of a pair of purely imaginay complex conjugate eigenvalues. Then lemma [I] 
and the fact that the Kronecker sum of two matrixes has all pairwise sums of eigen- 
values of the matrixes in the sum (see appendix |5| directly implies that the Jacobian 
matrix of the full system (Fni has the eigenvalues {±A, ±2A, 0} and the eigenvalue 



u>n = 15 W12 = —12 
w 2 i = 16 W22 = —5 
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is of multiplicity one. Therefore at this point, the Jacobian matrix of the infinite size 
system has its 5 eigenvalues having a null real part, and thus at this point the system 
is very degenerate and in particular the bifurcation is not a generic Hopf bifurcation. It 
is quite difficult numerically and analytically to study the solutions emerging from this 
bifurcation point, and many solutions can appear at this point, including cycles and 
fixed points. Similarly, at the two saddle-node bifurcations of WC system, since the 
Jacobian matrix of the infinite-size system has all the pairwise sums of the eigenvalues, 
it will have two null eigenvalues and therefore a degenerate bifurcation. The behavior 
of the system around these degenerate points will appear more clearly in the finite-size 
unfolding of these degenerate bifurcation points, and we will observe in particular that 
the Hopf bifurcation point corresponds to the merging of two generic Hopf bifurcations 
and the saddle node bifurcations to the merging of two generic saddle-node bifurcations, 
see section |2.2.2| and figure Fig. [6] In this section we do not address here the problem 
of classifying all behaviors of the system around the bifurcation (since as already men- 
tioned, the finite-size unfolding will help answering these questions), and restrict the 
analysis to the question of whether the WC cycle keeps a nonlinear stability in the new 
system. To answer this question, we numerically perturb WC's cycle by adding a small 
initial condition in the correlations. Interestingly, we observe that a new cycle appears, 
that is attractive in a certain range of parameters (see Figure [T]). The WC cycle does 
not completely loses stability and still appears as the only behavior of the system in a 
certain range of input values (e.g. I\ = —0.5 as plotted in Figure[TJb)). But for smaller 
values of the parameters, it loses its nonlinear stability and the additional cycle appears 
(e.g. h = -2 and Figuregc) and (d)). This new cycle is totally different in its shape, 
has a period close to half the period of the cycle corresponding to WC system, and has 
non-zero correlations that vary periodically at the same frequency. A continuation of 
this branch of cycles shows that it is stable in a significant range of parameters (Fig- 
ure [TJe)). It loses stability when I2 decreases through a period-doubling bifurcation, 
and as I2 increases through a Neimark-Sacker bifurcation. This branch of limit cycles 
emerges from the very degenerate point corresponding to the Hopf bifurcation of WC 
system, as we expected. It disappears at a point corresponding to a subcritical Hopf 
bifurcation with correlations (on a branch of unstable fixed points that is not plotted 
in the diagram but that will be further investigated in section 3.2 1. 

We conclude with a further observation on transient behavior. We have seen that 
WC's cycle was neutrally stable, but in some regions of parameters, it has nonlinear 
stability. In the region of parameters where the WC cycle is nonlinearly stable (white 
region in figure[TJd)), the transient phase of convergence towards this cycle takes a long 
time. In the green region where it is nonlinearly unstable, the neutral stability of the 
cycle produce some very strange transient behaviors where the activity of the cortical 
column is trapped around the neutrally stable WC system for long times, before leaving 
the neighborhood of this solution and converging towards the only stable solution which 
is the new cycle. 



Correlation-induced oscillations Correlations can have even more dramatic effects than 
modifying a periodic orbit as shown in Model I, and in this section we describe a case 
where a periodic orbit arises for parameter values where Wilson and Cowan system 
only fixed-points. To this end, we choose a two-population network known to have 
bistability. This model consists of two-populations with inhibitory interactions and no 
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(a) Bifurcation Diagram, (b) Wilson and Cowan cycle (c) New cycle of the infinite 
Wilson and Cowan system 




Fig. 1 Model I: Comparison between WC system and the infinite-size system: (a) and (d) are 
the bifurcation diagram with respect to the parameter Ii of respectively WC and the infinite- 
size system. The thick lines represent stable fixed point, thin lines unstable fixed points, plain 
circles the extremal values of a stable cycle and empty circle unstable cycles, (b) represents 
WC cycle as a solution of the infinite-size system, which shows null correlations, red: v\, white: 
i/2, yellow: An, green: A22, blue: A12 = Ai\. For non-zero correlations of the initial condition, 
the solution converges towards a new cycle with correlactions (c) and (d) which is stable in a 
determined region of parameters. 



self-interaction. The connectivity matrix in that case reads: 

/ w\\ = W12 = —12 \ 
\w 2 i = — 12 w 2 2 = ) ' 

The inactivation constants are assumed constant equal to 1, and the inputs to the 
two populations are distinct and denoted I\ and I2 (hence the system is not fully 
symmetrical). The firing rate function are also identical and we make the usual choice 
fi(x) — 1/(1 + exp(— a;)) for i £ {1, 2}. These parameters define a model called Model 
II. 

We break the symmetries between the two populations by considering that each 
population can receive different input, we freeze the input I2 and let I\ vary. In that 
case, the standard WC system presents a bistable behavior between two stable fixed 
points as we illustrate in Figure [2] The infinite size system driven by equations |7| 
with these parameters has much richer dynamics, as illustrated in Figure [2] We first 
set I2 = and study the bifurcation diagram as I\ varies. We observe that the infinite- 
size system features a stable limit cycle, which is a qualitatively nontrivial effect of 
the correlations in the mean-field limit. All fixed-points behaviors in WC system are 
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preserved, as predicted by theorem[T] (black curves in Figure[2jb)), but new fixed points 
with non-zero correlations appear (blue lines in the diagram). Along one of the newly 
generated branch of fixed point, a supercritical Hopf bifurcation appears, generating 
a family of stable limit cycles (red circles in Fig. |5|b)) that further undergoes few 
consecutive period doubling bifurcations. We note that the new branch of fixed points 
intersects the branch of fixed points of Wilson and Cowan system precisely at the 
saddle-node bifurcations of Wilson and Cowan system. This fact is compatible with a 
partial result of the proof of theorem [T] at the saddle-node bifurcation, instead of a 
single eigenvalue equalling zero, two eigenvalues are simultaneously equal to zero, and 
the saddle-node bifurcations of WC system appear as transcritical bifurcations in the 
infinite size system. Note that in that case, the WC behavior is nowhere the unique 
stable behavior of the system, and for any value of the input parameter I\ , the infinite 
system will present attractive behaviors that are different from the WC behaviors. In 
particular, it has stable fixed points with non-zero correlations and non-zero correlation 
oscillations. However, we note that though these behaviors exist in the dynamics, the 
newly discovered fixed points do not constitute acceptable solutions since they are 
characterized by a negative firing-rate V2 and a non-positive definite correlation matrix. 
However, the cycles constitute acceptable solutions that might have counterparts in the 
Markov system. 




Fig. 2 Model II, Wilson and Cowan (left column) and infinite-size (right column) systems, 
with I2 = (top row) and I2 = 5 (bottom row). 
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We further explore this model by choosing a large value of the input parameter 
to the second population, I2 ~ 5. For this level of input, the bifurcation diagram of 
WC system presents exactly the same features as in the case I2 = 0, namely two 
saddle-node bifurcations and a region of bistability of fixed points. The infinite-size 
system, similarly to the case I2 = 0, also has additional branches of fixed points with 
non-zero correlations that are connected to WC bifurcation diagram at the saddle-node 
bifurcation points. These saddle-nodes are degenerate and non-generic in the infinite 
size system from lemma[T]and become transcritical bifurcations. These additional non- 
zero correlation fixed points can be stable or unstable, depending on the parameter, and 
add additional possible behaviors. From one of these branches, the system undergoes a 
supercritical Hopf bifurcation and a family of stable limit cycles appears. This branch 
of limit cycles undergoes two period-doubling bifurcations, torus bifurcations and two 
fold of limit cycles, yielding additional oscillatory behaviors to the purely fixed-point 
structure of WC system. 

From all these observations on both models, we conclude that the infinite-size 
system, besides featuring the same solutions as Wilson and Cowan system as a subset 
of solutions with zero correlations, also has additional behaviors that qualitatively 
differ from the Wilson and Cowan system. The mean-field limit including correlations 
is therefore a more complex system than WC system, in which correlations of the 
firing activity interact in a non-trivial way with the mean firing rate yielding complex 
behaviors. These qualitative correlation-induced behaviors can therefore be a good 
indicator of which model better accounts for the behavior of cortical columns and 
cortical areas if these behaviors hold in the finite-size systems which are equivalent to 
the moment equations of the initial Markov chain. 

We now turn to investigating finite-size effects arising from the finiteness of the 
number of neurons in a network. 



3 Finite-size effects 

The bifurcation diagram displayed by the infinite-size system studied in the previous 
section is composed of some very degenerate points, and displays behaviors the WC sys- 
tem does not display. The question that naturally arises at this point is to understand 
whether the infinite-size model or WC model really capture the qualitative behavior of 
large networks in the Markovian framework through the BCC and Bressloff finite-size 
equations. 

In this section we first study how the infinite-size bifurcation diagram is unfolded 
in finite-size networks, before simulating the Markov process and comparing its global 
behaviors to the behaviors of the finite-size systems. We observed that the finite-size 
BCC and Bressloff equations appeared as a perturbation of the infinite-size system, in 
which the size of the network appears as a parameter of the equations. In all of this 
study, instead of considering discrete population sizes, we will consider the inverse of 
the total population size n = 1/N as a continuous parameter. Non-trivial distinctions 
between finite-size networks and infinite size equation will essentially appear at the 
bifurcation points. Indeed, since both the rescaled Bressloff and BCC models vector 
fields are continuous with respect to the parameter n, when the infinite-size system 
(n = 0) exhibits an hyperbolic fixed point, this fixed point will have a counterpart with 
the same stability for finite, sufficiently large networks, by application of the implicit 
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functions theorem. This is also the case of hyperbolic fixed points of WC system, since 
these are hyperbolic fixed points of the infinite-size system with zero correlations. In 
order to study similarities and qualitative differences between finite-size networks and 
the infinite size limit, we will therefore focus on the bifurcations of both systems and 
numerically study the codimension two bifurcations in the one and two populations 
cases, with respect to the size of the network and to another parameter of interest (for 
instance an input parameter). 



3.1 Analysis of single-population finite-size networks 

We start by investigating a single-population case. In the mean-field limit, we have 
shown that the infinite size system essentially exhibits the same behaviors as the WC 
model in terms of fixed points and limit cycles. But is this still the case where the 
number of neurons is finite? 

At the bifurcations of the WC system, singular phenomena can occur. Consider for 
instance a saddle-node bifurcation point of WC system. At this point, the infinite size 
system will present a double zero eigenvalue, and the unfolding of this bifurcation can 
lead to different cases. As we show in this section, it can either unfold into a generic 
saddle-node bifurcation as we show in the first section, or it can generate an oscillatory 
behavior through a limit cycle, creating in the finite-size system oscillations which do 
not exist in the WC or in the infinite size one population equations. 



3.1.1 Regular unfolding of the infinite-size system 

We start by considering BCC model with a strictly positive voltage to rate function 
/, and focus on acceptable solutions that have positive firing rates and correlations. 
In that case, we numerically show that there is no qualitative difference between the 
stable solutions of the finite-size BCC network on one hand, and the infinite size system 
and WC system on the other hand (we showed that both models present the same 
qualitative behaviors in the one population case). We address this problem numerically, 
with an inactivation constant a = 1 (without loss of generality, since modifying this 
value only amounts rescaling time), w — 10, and letting both the size of the network 
n and the input parameter I free. The activity to rate function is chosen to be f(x) — 
l/(l+e*p(-a;)). 

With these parameters, WC system presents two saddle-node bifurcations when 



the input parameter / varies, as shown in Figure 3(a) The system therefore presents 
an interval of values of the parameter I where the system presents two stable fixed 
points and an unstable fixed point, and outside this interval, the system has a unique 
fixed point which is stable. Similarly to what we showed in the previous section, the 
infinite-size network presents exactly the same features, and no cycle. 

We now study BCC finite-size equations. We set the size of the network to be 
N — 50. BCC finite size equations present the same qualitative behaviors, and the 
same type of bifurcations: the system undergoes two saddle-node bifurcations as the 
input parameter value / varies, and a bounded interval of input values / where the 
system presents a bistable behavior. The only qualitative distinction between finite- 
size networks and infinite size networks is the existence of a two unstable fixed points 
in the bistable zone and the fact that the system presents two distinct branches of 
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fixed points, but these differences only affect unstable fixed points and do not produce 
generic qualitative differences in the behaviors. 

The way this bifurcation diagram is transformed into the bifurcation diagram of 
the infinite system is very simple. The two distinct branches of fixed point get increas- 
ingly close, and in the limit n = 0, the two branches are superimposed. For finite-size 
networks, we observe that the value of the correlations of the stable fixed points was 
non-zero (see Fig.[3|. These values were nevertheless quite small, and decrease towards 
zero as the size of the network as the number of neurons increases. When continuing 
the two identified saddle-node bifurcations as the network size increases, we observe 
that the bifurcation curve persists until the limit n — with no singularity. At the limit 
however, the non-zero eigenvalue of the Jacobian matrix at the saddle-node bifurcation 
point tangentially reaches the value yielding the singularity predicted in proposition 
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(c) BCC bifurcations (Correlations) (d) Codimcnsion two bifurcations 



Fig. 3 BCC model and the infinite-size system (a) Bifurcation diagram of WC system with 
respect to the input /. Plain curve: stable fixed point, dotted curve: unstable fixed points, 
LP stands for limit point (or saddle-node bifurcations) (b) and (c) Bifurcation diagram of 
BCC model with respect to I, for N = 50 neurons, (b) Firing rates, (c) Correlations, (d) 
Codimension two bifurcation diagram with respect to I and the network size n = 1/N. Pink 
line: Saddle- node bifurcations (Limit Points). The star at n = denote double zero eigenvalue 
in the Jacobian matrix singular point. The values correspond to the values of WC saddle-node 
bifurcations. 



In this case, we therefore illustrated the fact that BCC system can smoothly unfold 
the WC system by introducing correlations. These correlations vanish as the size of 
the population increases, yielding a purely Poisson firing with firing rates given by the 
solution of WC equations. For the same parameters and functions, it is easy to show, 
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following the same steps, that Bressloff rescaled model also smoothly unfolds the WC 
system and that the finite-size system presents exactly the same type of bifurcations 
with two separated branches of fixed points each of which collapse in the limit where 
the number of neurons is infinite (not shown), and the value of the correlations of the 
fixed points decreases towards zero when the number of neurons increases. However in 
this case, one has to bear on mind that the interpretation of the result slightly differs: 
the limit behavior is now seen as a purely asynchronous state where all the neurons 
of the population fire independently in the limit TV — > oo. These behaviors will closely 
match the Markovian model's behavior as shown in section [3. 1 . 1[ In these cases, WC 
system conveniently summarizes the behavior of large networks, and provides a simple 
model to study large networks. However, there exist models for which this conclusion 
does not hold, and where the infinite-size system's behavior unfolds into more complex 
behaviors. 

3.1.2 Singular unfolding of the infinite-size system 

We now turn to study the possible singular unfolding of the infinite size system in one 
population. To this purpose, we will be specifically interested in studying oscillatory 
behaviors in the finite-size systems that are not present in the limit where the number of 
neurons tends to infinity. We will call these qualitative distinctions between finite-size 
and infinite-size systems singular unfolding. 

To this purpose, we choose a framework where fixed points are easily described, 
and along the fixed point search for sufficient conditions to get a supercritical Hopf 
bifurcation, yielding oscillatory behaviors in single population networks that are neither 
present in the infinite-size system nor in the WC system. To fix ideas, similarly to the 
previous case treated, we consider BCC one population model. We consider here a 
non-physiological cases by allowing the firing rate variable to be negative. We choose 
an activation function / such that /(/) = for some fixed J, which can always be done 
with an unspecified sigmoid / through the modification f(x) = f(x) — /(/). This quite 
unrealistic case will allow studying the qualitative differences that can appear between 
the finite and infinite size systems. However, the fixed point is quite particular, and 
we will see that this phenomenon is not generalized in a simple way when considering 
positive activation functions / and positive firing rates (the affine transform f(x) — f(I) 
does not have a trivial effect on the nonlinear system). 

The BCC single population model reads: 

^ = -av + f(wv + I) + \f"(wv + I)w 2 C 

Therefore, under the assumptions made, it is easy to see that the null solution v — 
0, C — is always solution of the system. The Jacobian matrix J at this point reads: 

f-a + wf'(s) \f"{s)w 2 \ 
\ 2/'( S )$ 2 {-a + wf'(s))J' 

where we denoted s = w u + I — I for v = 0. The trace of the Jacobian matrix at this 
point is equal to 3(— a + w f'(s)) and the determinant is equal to 2(— a + w /'(«)) — 
f'(s)f"(s)w 3 /N. The trace vanishes at the null fixed point for (a,w) such that 



a = wf'(I), 
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and therefore there exist a set of parameters for which the trace vanishes. In order to 
ensure that this null trace corresponds to the presence of a pair of purely imaginary 
eigenvalues, we further need to ensure that the determinant is strictly positive. When 
the trace vanishes, the determinant simply reads —f'(s)f"(s)w 3 /N. Since / is a sig- 
moidal function, /' is always positive, and therefore the determinant is positive as soon 
as / (s) < 0, which can be ensured for some specific firing rate functions. In that case, 
we will denote ujq = \J — f (s) f" (s)w 3 /N . 

Therefore, under the assumptions that there exists / such that /(/) = and 
f"(I) < 0, there exist parameters a and w such that the Jacobian matrix of the 
system at the null solution has a pair of purely imaginary eigenvalue. Let us now check 
that this point corresponds to a Hopf bifurcation when varying the parameter a. To 
this purpose, we check transversality and genericity conditions of the Hopf bifurcation 
at this point, as given for instance in |35| . These calculations are provided in appendix 
[C| The transversality condition is easily checked. It consists in showing that the dif- 
ferential real part of the eigenvalues of the Jacobian matrix at this point, with respect 
to the parameter a, does not vanishes. Let us denote by fi{a) the real part of the 
eigenvalues. Since we are in a planar system, we have: 

/*(") = \tt{J) 

dM(oQ = _3 < Q 
da 2 

Therefore the transversality condition of the Hopf bifurcation is satisfied at this point. 
In order to fully state that the system undergoes a Hopf bifurcation, we need to show 
that the first Lyapunov exponent of the system at this point does not vanish. This 
Lyapunov exponent, after some tedious calculations using the formula given in [35] 
(see appendix [c|) , is shown to be of the same sign as: 

>>(j)/(D (i + ^) +/ "(/) 2 (^^Y . 

The sign of this expression is governed by the values of the differentials of / at the point 
I. If the expression in the bracket is negative, then the Hopf bifurcation is supercritical 
and the system has stable oscillations. Let us now investigate the dependence of this 
cycle on the network size. We have shown that ujq = \J — f (s)f' (s)w 3 /N . Therefore, 
if this cycle appears for some finite network, it will lose stability as the size of the 
population increases. The period of the generated cycle is proportional to the inverse 
of the square root of N, and therefore slowly tends to zero as the network size increases 
(note that the cycle will lose stability as the network size increase, since for sufficiently 
large populations, the Hopf bifurcation will be subcritical). 

Now that we identified simple conditions for the Hopf bifurcation to take place, we 
exhibit a network that indeed has this bifurcation. We choose for instance as activation 
function an hyperbolic tangent function 



w z N 
JW 



f(x) = tanh(a;) — tanh(J) 
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with / = 0.5. In that case, we have: 



(f(I) 
/'(/) 
/"CO 

lf"'(I) 



= 

« 0.786 

» -0.727 < 

« -0.56 



Therefore this function satisfies the assumptions done in the previous paragraph. More- 
over, since /"'(/) < and (2/cjo - ^) < ( for TV small enough), the first Lyapunov 
coefficient is negative, and therefore the system undergoes a supercitical Hopf bifurca- 
tion and presents oscillations, as illustrated in Figure [4] 




(a) Oscillations 




(b) Continuation and Homotopy 



Fig. 4 Oscillations in a single population network with finite size, (a) is a cycle observe for /o 
(p=0), (b) is the codimension two continuation with respect to w and p of the Hopf bifurcation 
for different values of a: pink for a = .3, blue: a = .8 and green: a = 3. 



Therefore, we exhibited a case where the finite-size equations have important qual- 
itative differences from the infinite-size equations. However this system involved non- 
physiological values of the parameters, in particular a non-positive activation function 
/. Therefore, this behavior will never appear in the Markov process that yielded BCC 
equations. 

The question that then arises is whether such behaviors appear in plausible set- 
tings. In order to answer this question and search for a system with positive sigmoidal 
functions that indeed present oscillations in one-population networks, we continued the 
Hopf bifurcation while continuously transforming the non-positive sigmoidal transform 
into a positive sigmoid. For instance, we consider / a non-positive sigmoidal function 
for which the Hopf bifurcation exhibited here appears, and define f p = f — pmin(/) 
a transformation of / such that /o = / and /i is a positive sigmoid. This homotopic 
transform smoothly maps the vector field presenting an attractive cycle onto a vec- 
tor field with a positive activation function. Since we have proved the genericity and 
transversality conditions on the Hopf bifurcation identified, we are in position to con- 
tinue it as the parameter p varies together with a in a codimension two bifurcation 
diagram (by application of the inverse function theorem). We observe in that case that 
the Hopf bifurcation cannot be continued when varying the homotopy parameter until 
reaching a positive sigmoid (see Figure Rib) for the example provided above). We have 
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tried numerous other sets of parameters and have never been able to continue periodic 
orbits to a regime where there is a non-negative activation function. 

Therefore, though the one-population BCC equations can display oscillations for 
non biologically relevant parameters, we conjecture that there is no acceptable oscilla- 
tory solution for positive activation functions. 

3.2 A two-populations finite-size network 

We have seen in the previous section that the one-population finite-size BCC model was 
accurately approximated by the infinite-size system (and therefore the WC system from 
the result of a previous section) for biologically realistic parameter, though oscillatory 
phenomena may occur for non-biologically relevant parameters that allow both the 
firing-rate and the correlations to be negative. The picture will be quite different in 
multi-population networks. Indeed, in higher dimensions, the WC system can have Hopf 
bifurcations that create very degenerate bifurcation points in the infinite-size system, 
which may non-trivially unfold in the finite-size systems. In this section, we consider 
the two population networks with excitatory self-interaction and negative feedback 
loop, called Model I, introduced in sect ion [2 . 2 . 2 1 and address first the case of Bressloff 
model, before showing that BCC model presents the same features in appendix [D] 
This study will allow addressing in particular the question of the existence of so-called 
quasicycles and also of aperiodic solutions in the system close to the Hopf bifurcation 
of the mean-field Wilson and Cowan limit, addressed recently in the same model by 
Bressloff in These quasicycles correspond to intrinsic noise-induced oscillations 
in parameter regimes where the deterministic mean-field limit (in that case Wilson 
and Cowan system) is below the Hopf bifurcation, and where the system features 
a stable focus. This phenomenon was also recently found in similar systems arising 
in biochemical oscillations of cellular systems, see e.g. [9ll39]. Our approach to these 
phenomena consists in unfolding the Hopf bifurcation of the mean-field system to 
finite network sizes, in order to understand how the intrinsic noise (that can be seen 
as through the Langevin approximation of the system valid for large population sizes 
as used in appendix [A] and in [11]) interacts with the Hopf bifurcation. 

To fix ideas, we focus in this section on the particular network introduced in sec- 
tion |2.2.2| namely the two-population rescaled Bressloff model. Denoting c = C, the 
equations read: 

' v[ = -aivi + f(si) + 5 /"(si) (wfi en + wj 2 c 22 + 2 iu 12 Wu c i2 ) 

v' 2 = -OL 2 V2 + f(s2) + \ .f"{s2) (»22 c 22 + U>21 C H + 2 ™22 ™21 C12) 

c 'n = w[[ a i u i + /Oi)] - 2aicn + 2/'(si) (wu en + w 12 c.12) 

c 22 = TfJ [ Q 2^2 + /(S2)] - 2a 2 C22 + 2/'(s 2 ) (w 2 l C12 + W 22 C22) 
c 12 = -( a l + Q 2)C12 + /'(«l) (mi C12 + W12 C22) 
+f(s2) («J21 Cll + W22 C12) 

where we denoted: 

J s\ = w\\ v\ + W12 v % + h 

|^S2 = W21 v 2 + W22 V2 + *2 

For the set of parameters studied in this section and a number of neurons iV = 50, 
Bressloff finite-size system undergoes four Hopf bifurcations labelled HI through H4 
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and six saddle node bifurcations LP1 through LP6. It features two distinct families of 
limit cycles, one that undergoes two folds of limit cycles LPC1 and LPC2, and one 
that undergoes a Neimark-Sacker (Torus) bifurcation, as displayed in Figure [5^b), and 
are compared with the standard WC system that displays one Hopf bifurcation point 
H and two saddle node bifurcations LPa and LPb. 




(a) Wilson and Cowan system (b) Bressloff model 



Fig. 5 Bifurcation diagram of Model I: (a) Wilson and Cowan system and (b) Bressloff sys- 
tem. Blue: equilibria, pink: cycles. Bifurcations of equilibria are denoted with a red star, LP 
represents a saddle- node bifurcation (Limit Point), H a Hopf bifurcation. The four Hopf bifur- 
cations share two families of limit cycles. One of these undergoes two fold of limit cycles LPC, 
and the other branch of limit cycle a Neimark Sacker (Torus) bifurcation. 

The analysis of this bifurcation illustrates the fact that the small finite-size system 
studied have clear qualitative differences with the WC system. However, because of the 
properties of the vector field, and in particular its smooth dependency with respect to 
the size of the network N, the bifurcation diagram of Fig. [5] smoothly connects to the 
degenerate bifurcation diagram of the infinite-size system presented in section [2. 2. 2| 

Let us now investigate how the bifurcation diagram of the infinite system unfolds 
into the bifurcation diagram of the finite-size system. To this end, we continue of the 
bifurcation diagram of figure Fig. [5|b) as the number of neurons increases, as shown 
in the codimension two bifurcation diagram of figure Fig. [6] The structure of the new 
bifurcation diagram appear way more intricate than WC's, and in particular, each 
bifurcation point of WC system, corresponding to degenerate dynamics of the infinite- 
size system, unfolds non-trivially into different generic bifurcations in the finite-size 
moment equations. Let us describe the diagram plotted in figure Fig. [6] We observe 
that the bifurcation diagram of the finite-size system, for any N < 232 neurons will 
be similar to the one plotted in figure Fig. [5] As TV increases, the bifurcation diagram 
is smoothly transformed except at two particular codimension two bifurcation point: a 
cusp bifurcation at n « 0.0867, from which emerge two additional saddle-node bifur- 
cations denoted LP7 and LP8. The value of the correlation variables of these points 
progressively converge to zero as the network size increases. These two saddle-node 
bifurcations correspond in the infinite-size limit to the two saddle nodes of WC system 
LPa and LPb, and are not present in the finite-size model for a small number of neu- 
rons. The second codimension two bifurcation appearing in the diagram correspond 
to the disappearance of the Hopf bifurcation point H4 through a Bogdanov-Takens 
bifurcation with the saddle-node bifurcation manifold corresponding to LP8. Except 
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from these two codimension two bifurcations, the continuation of all other bifurcations 
is smooth for any n > 0. But in the limiting case n = 0, the system is highly singular 
and different bifurcation manifold collide, as we now discuss in more details. Let us 




Fig. 6 Codimension two bifurcations with respect to the input to population 1 and the popu- 
lation size in Bressloff 's model. Pink: Hopf bifurcations, blue: saddle node bifurcations, green: 
folds of limit cycles. The continuation of the Hopf bifurcations HI and H4 merge at WC's H 
point when N — > oo. We observe that the fold of limit cycles merge with other bifurcations in 
the limit of infinitely many neurons (see text). In the limit N — > oo (corresponding to n = 0) 
we labelled the points H, LPa and LPb arising in WC system. Points with zero correlations 
are circled in red, all other points have non-zero correlations. 

first follow the continuation the Hopf bifurcation labelled HI in diagram of Fig. [5|b) 
as the network size increases. We observe that this Hopf bifurcation collides, in the 
limit N — > oo, with the continuation of the Hopf bifurcation H3 of the second branch 
of fixed point. Their collapse corresponds to the degenerate Hopf bifurcation point of 
the infinite model, denoted H in Fig.[5|a). The cycle originating from HI and the cycle 
originating from H3 will therefore emerge from this fixed same point, accounting for 
the result presented in section [2] and in particular figure Fig. [I] where we evidenced the 
presence of two different cycles emerging from the same bifurcation point H. Moreover, 
we have seen that the corresponding branch of limit cycles in the infinite-size system 
bifurcation diagram Fig. [l] loses stability. This loss of stability corresponds to the con- 
tinuation of the Neimark-Sacker bifurcation NS of figure [5] that can be continued in 
the mean-field limit. However, an important distinction between the finite and infinite 
case is that This cycle undergoes an additional period-doubling bifurcation when the 
size of the network increases, as observed in the bifurcation diagram of the infinite-size 
system, and as shown in figure [7] In that figure, we continued these cycles for one fixed 
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value of the parameter I\ = —0.5, as the number of neuron increases. The family of 
limit cycles continued is stable in a bounded interval of values for the parameter N. 
For a number of neurons close to TV = 142 (precisely ri2 = 0.01418), the branch of limit 
cycles undergoes a fold of limit cycle, and therefore the branch of limit cycles does not 
exists for smaller networks. For networks larger than N = 142 neurons, the system 
presents a branch of stable and a branch of unstable limit cycles. Stable cycles persist 
up to networks as large as JV = 1700 neurons (precisely n2 = 0.001174). At this size, 
where the branch stable limit cycles changes stability through a subcritical Neimark 
Saker (Torus) bifurcation with no strong resonance. When increasing further the size of 
the network, the system shows the presence of a chaotic attractor which clearly cannot 
exist in the two-dimensional WC system, yielding an important qualitative distinction 
between large scale networks and the relative mean-field limit (see [7| . This chaotic 
behavior exists until the infinite-size is reached. 




Fig. 7 Chaotic behavior for large networks: (a) bifurcation diagram as a function of the 
network size. NS: Neimark-Sacker (Torus) bifurcation and FLC: fold of limit cycle, (b) Poincare 
maps representing ai and en with Poincare section a,2 = 0.7. 



Let us now follow the branch of limit cycles arising from H3 as the network size 
varies. We observe in our codimension three bifurcation diagram that the branch of 
Hopf bifurcations H3 exists until the infinite-size is reached. In that limit, we observe 
that this branch of Hopf bifurcation asymptotically meets the fold of limit cycles LPC2 
and the curve of saddle-node LP7 and the saddle homoclinic curve arising from the 
Bogdanov-Takens bifurcations. This point corresponds therefore to a very degenerate 
behavior of the system, and corresponds, in WC bifurcation diagram, to the saddle-node 
bifurcation LPa that also corresponds to the homoclinic bifurcation point. Regarding 
the fourth Hopf bifurcation H4, we already observed that it was continued until it 
disappeared, in finite size, through a Bogdanov-Takens bifurcation. 

Let us now consider the continuation of the saddle-node bifurcations. The bifurca- 
tion diagram of WC presents two saddle-node bifurcations LPa and LPb that corre- 
spond to the two saddle-node arising from the cusp bifurcation CP of Fig[6] that are 
labelled LP7 and LP8. In the mean-field limit, the saddle-node bifurcations correspond- 
ing to the continuation of LP1 and LP4 merge in a cusp bifurcation and disappear. At 
this same point, the fold of limit cycles LPC1 collapses and disappears in a homoclinic 
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bifurcation. The saddle-node bifurcations LP2 and LP6 disappear in the mean-field 
limit by being continued tangentially to the line n = 0. The continuations of the saddle- 
node bifurcations LP3 and LP6 crosses the line n = with no singularity and therefore 
exists in the infinite size system, corresponding to a non-zero correlation saddle-node 
bifurcation and that do not correspond to any plausible behavior (all plausible be- 
haviors correspond to zero correlation points). Eventually, the saddle-node bifurcation 
LP7 connects with the fold of limit cycle, generating the homoclinic bifurcation of the 
infinite-size system. 

We therefore observe that the bifurcation diagram of the finite-size Bressloff model 
non-trivially unfolds the infinite-size system, and that at any finite scale the bifurcation 
diagram and the solutions it present differ significantly from the infinite-size system 
and from WC system. In particular, it shows additional cycles, as observed in the 
infinite-size system previously, and a chaotic effect that is not present in the two- 
dimensional WC system. This analysis also illustrates a mesoscopic-scale phenomenon 
in the presence of a periodic orbit that is stable only when the size of the network is 
comprised roughly between one hunded and one thousand neurons, a typical size of a 
cortical columns often modelled by Wilson and Cowan system. These results remain 
valid for BCC model, as shown in appendix[D] The question that now arises is whether 
these behaviors indeed reflect behaviors of the Markovian model. 



3.3 Finite-size Markovian model 

In the previous section, we investigated the distinctions between the solutions of the 
finite-size BCC model and compared these with those of the infinite-size and of WC sys- 
tems. We now return back to the original Markovian system that allowed derivation of 
Bressloff and BCC equations, and compare simulations of this network with the finite- 
size Bressloff model, composed of two differential equations arising from the moment 
truncation of this random variable. We will perform this analysis in the two models 
with two populations investigated in the paper, Model I and Model II. This analysis is 
particularly interesting in regions of the parameter space where the system presents a 
fixed-point behavior, and where the moment expansion is a relevant representation of 
the system. In the periodic orbit regimes, we discuss a suitable methodology, beyond 
the scope of the present paper, that is suitable to analyze the system. 

In order to simulate the Markov chain, we make use of Doob-Gillespie's algorithm 
|20II30||3T] . From the master equation governing, e.g. BCC model: 

M 

— dt = [ a * ( ni + i ) p ( n '+> *) - Q * n * p ( n > *) 
i=i 

+ Fj (nj_ )P(m- ,t)-Fi (n) P(n, t)] 

we derive the transition rates of the state variable n. We use these rates to simulate 
sample paths of the state variable n and rescale these to compare to the solutions 
obtained for BCC or Bressloff models and the infinite-size models. Given the con- 
figuration n(t) of the network at time t, we have in the BCC case (Bressloff case is 
straightforwardly treated using the same method): 

— The probabilistic intensity for the transition rii(t + 1) = rij(i) — 1 is equal to 

Qi =airii, 
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— The probabilistic intensity for the transition m(t + 1) = n;(t) + 1 is equajj to 

Given an initial condition n(0), a simulation consists of computing all possible tran- 
sitions probabilities and fi. We then draw the time for the next transition as an 
exponential random variable with intensity Q = ^ife + fi)- This time provides the 
next event occurring in the chain. By the properties of exponential laws, the transition 
is rii(t + 1) = rii(t) — 1 (resp. rn(t + 1) = rn(t) + 1) for some i £ {1, . . . , M} with 
probability qi/Q (resp. fi/Q). This simulation algorithm for the Markov process n(t) 
is exact in law, and does not involve any approximation such as time-discretization, 
and therefore allows efficient simulation of the sample path and of its statistics. These 
statistics are computed using a Monte-Carlo algorithm consisting of simulating many 
independent sample paths of the process and deriving from these the trajectory of 
the mean firing rate i>i(t) and the correlation in the firing activity Cij(t) through the 
empirical mean and correlations obtained. In the simulations we focus on BCC model. 
The same numerical experiments can be performed in the Bressloff case. 

3.3.1 One-Population Model 

We start by considering the one-population model studied in section [3. In that case 
we have shown that both Wilson and Cowan and the infinite-size systems had the same 
behavior, and that this behavior was regularly unfolded in the finite-size systems. In 
these cases, two saddle-node bifurcations defined a zone of bistability that depended on 
the size of the network as shown in Figure |3^d) . We simulate the network using Doob- 
Gillespie algorithm and identify the boundaries of the zones where bistability occurs 
by continuation of the fixed points obtained as the parameters increase or decrease. 
The continuation is initialized by running the algorithm with parameters associated to 
a single attractive fixed point of WC system. Then the parameter / is slowly varied 
and simulations are run with initial condition being the last values of the network state 
of the previous simulation. We observe the hysteresis phenomenon and bistability, and 
the bistability zones closely matches WC system, as shown in Figure|8] We simulate the 
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Fig. 8 Bistability zones in BCC finite-size model (black dashed curve) and computed from 
the Markovian BCC model (colored plain curve). 

network for increasingly large networks ranging from N = 50 neurons to N = 100 000 
neurons with 10 000 simulated sample path from which we extract the empirical mean 




1 Note that the activation functions are different in BCC and Bressloff models. 
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and standard deviations. From the curves obtained in that fashion, we automatically 
identify zones where there is bistability. These boundaries are plotted together with 
the values of the saddle-node bifurcation at these points arising from the codimension 
2 bifurcation diagram of BCC finite-size model, and we observe that the curves closely 
match as soon as the network is larger than 1 000 neurons. Therefore in this model, BCC 
system closely reproduces the behavior of the finite-size Markov model for sufficiently 
large networks. A similar result is obtained in BCC case (not shown). 

3.3.2 Markovian Model I : Fixed points regimes 

We now turn to study the two-dimensional network Model I. In that case we showed 
that depending on the input parameter I\ , Wilson and Cowan system can either present 
a stable cycle or a stable fixed point. We first chose a case where WC system converges 
towards a fixed point, and to this end will study the system for I\ = —5. In that case, 
both WC system, the infinite-size system, BCC and Bressloff systems converge towards 
a fixed point. The periodic regimes are addressed in section [3. 3. 4| 

In the case where WC system presents a unique hyperbolic attractive fixed point, 
we showed that the infinite-size presents an exponentially stable fixed point with the 
same values of mean firing rate and null correlations, which has a counterpart in the 
unfolding of BCC and Bressloff model. Simulations of the Markovian model illustrate 
the fact that all the sample paths of the process converge quite fast towards this fixed 
point with null correlations displaying stochastic variations, and each sample path 
presents a precise match with the solution of WC solution at this point. The stochastic 
variations around Wilson and Cowan trajectory arising from the randomness of the 
Markov chain are averaged, and produce a fixed mean firing-rate corresponding to the 
fixed point of the infinite-size system, i.e. having null correlations and the same value of 
mean firing rates. Moreover, we observe that the correlations vanish, which implies that 
the state obtained is an asynchronous state, as expected from Bressloff 's expansion. 

Figure [9] represents simulations of Model I with I\ = —5. In panels (a) and (b) are 
plotted four sample paths of the process for 10 000 neurons simulated over on a time 
interval equal to 100, and panels (c) and (d) the mean-firing rates averaged over 5 000 
realization of the process, and thats shows a very close agreement with the solution 
of WC system. The maximal absolute value of the correlations (not plotted) are of 
the order of 10 -5 and are neglected as an effect of the finiteness of the sample path 
simulated for computing the mean and covariance. The same observations hold for all 
the simulations performed for parameters corresponding to an hyperbolic fixed point of 
the infinite system. Note that the difference between mean-firing rates of the solution 
of WC system on one hand and of BCC system on the other hand is of the order of 
10 -5 at this network scale, and therefore both models fairly approximate the stochastic 
Markovian model. 

3.3.3 Markovian Model II 

We perform the same analysis in the case of the network Model II studied in section 
|2.2.2| with symmetric connectivity with null diagonal and inhibitory interaction in BCC 
case. In that case, we showed that WC system presented bistability, and the infinite- 
size system presented a more complex behavior with multistability and different limit 
cycles. In the simulations of the Markovian system, similarly to the case of Model I, we 
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Fig. 9 Network Model I, for parameters associated with a stable fixed point (7i = —5) closely 
matches WC trajectory with random variations around the fixed point that average at the 
value of the fixed point, and correlations vanish (asynchronous state). (a),(b): The different 
colors correspond to different sample paths. The simulations were performed in the BCC case. 



observe that when WC system presents a single hyperbolic fixed point, the Markovian 
system closely matches with the solution of WC system during the transient phase, 
and the stationary state randomly varies around the fixed point of WC system. The 
mean of the firing rates over many realization of the process converges exactly towards 
the solution of WC system, and the correlations tend to zero, implying the system is in 
an asynchronous state and that its behavior is efficiently summarized by Wilson and 
Cowan system. We note that in that case, the infinite-size system and the finite size 
systems both presented additional fixed points with non-zero correlations, that do not 
actually exist as a regime of the Markov chain. 

For parameters associated with bistability in WC, the situation is similar, and 
we observe in our simulations that the WC system captures precisely the behavior of 
the Markovian system and all additional solutions derived from the moment hierarchy 
truncation do not appear. Simulations of the Markovian system show that depend- 
ing on the initial condition, the system will converge towards one or the other fixed 
point. Albeit the stochasticity of the Markov chain yielding the possibility to switch 
between the two fixed points, this phenomenon does not occurs often, and except for 
input parameters very close to the saddle-node bifurcations, no switching between fixed 
points is observed. In that case, we note that the additional fixed point with non-zero 
correlations and the additional cycles evidenced in the infinite-size system and that 
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are present in the finite-size moment equations do not appear in the simulations of the 
underlying Markovian model. We present in Figure [lO]the values of taken by the mean 
firing rate and correlations for 1000 realizations of the Markov chain over the 20 last 
units of time for a simulation of 20 000 neurons as a function of the input paramter I\ . 
We treat the cases I2 — and I2 ~ 5, the blue crosses are the continuation of the fixed 
points with I\ increasing (from —10 to 10 for I2 = and from —5 to 15 for I2 = 5) and 
red circles for I\ decreasing. The continuation is initialized by running the algorithm 
with parameters associated to a single attractive fixed point of WC system. Then the 
parameter I\ is slowly varied and simulations are run with initial condition being the 
last values of the network state of the previous simulation. We observe the hysteresis 
phenomenon and bistability, and the bifurcation diagram obtained perfectly matches 
WC's, and no counterpart of any of the identified correlations effects are observed. 
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Fig. 10 Stationary values of the Markov chain in Model II with I2 = as a function of Ii 
(see text). Continuation of equilibria by increasing I\ (blue crosses) and by decreasing I\ (red 
circles), superimposed with WC bifurcation diagram (black curve). Correlations are plotted 
with a different axis origin in cyan (Cu), magenta (C22) and yellow (C12), keep very close 
to zero except at the precise value of the bifurcations, where the jump between the two fixed 
points produces a spurious standard deviation. BCC case. 



It appears from this analysis that the finite-size moment equations present ad- 
ditional solutions, associated with large correlations, that do not exist in the initial 
Markovian model. We emphasize the fact that though the infinite-size model is not 
equivalent to the moment equations derived by both authors, the solutions evidenced 
have a counterpart in the finite-size moment equations corresponding to both systems. 
These equations present solutions with large correlations (that actually diverge as the 
number of neurons tend to infinity), which are not solutions of the initial Markov chain, 
which naturally behaves as WC system for large population sizes, as expected from the 
moment equations by both authors. 

A complementary approach to the study of bistability in the Markov framework al- 
lowing to quantify the probabilistic intensity of noise-induced transitions can be driven 
through the use of WentzelKramersBrillouin (WKB) approximation (see e.g. [211111] 
and references therein) performed directly on the Markovian system. This method was 
used by Bressloff in [11] in the same context for a one-population model, and more 
generally allows studying noise-induced transitions between multiple metastable states 
where the mean-field theory based on the system-size or loop expansions break down. It 
allows reducing the problem to a Hamilton- Jacobi equation from which one can derive 
the escape rate, namely the probability intensity that the solution escapes from the 
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attraction bassin of a given equilibrium of the deterministic equation (the WC system) 
and enters the attraction bassin of the other equilibrium. 

3.3.4 Periodic regimes 

The models proposed in [16il7l[TQ] emphasize the role of the mean-firing rate and of 
the zero timelag correlation functions. As illustrated in the previous section, these 
variables are suitable to model the system at fixed points. However, it remains unclear 
how these variables characterize the behavior of the Markovian system in periodic or 
chaotic regimes. 

The analytical study of such regimes can be carried out in a satisfactory fashion 
either in the theory of master equations and in the theory of stochastic dynamical sys- 
tems, using the Langevin equation. Such studies have been addressed by various authors 
within the context of chemical master equations [39,9,29 and recently Bressloff applied 
the same tools to neural models, using Langevin equations and Floquet theory [llj . 
Indeed, it appears that in regimes where the Wilson and Cowan system present a cycle, 
the Markovian model is characterized by random fluctuations in the neighborhood of 
this attractor, and the system loses phase information over time. These dynamics can 
be suitably studied in terms of a diffusive phase variable. In details, under the Langevin 
approximation of the dynamics, a possible way to study the behavior of the solutions 
of the Markovian system would be to use the stochastic phase reduction (see |48| ) . 

In this section we do not make use of these more customary methods but rather 
investigate the relationship between behaviors of the finite-size moment equations and 
behaviors of the Markovian system. For our choice of parameters in Model I, let us 
consider the parameter region I\ £ [I\ Hi II LPa] where I\ h (resp. I\ LPa) correspond 
to the value of I\ at the point H (resp. LPa) . In this region, the solutions of WC system 
are periodic orbits. The question we address is whether the behavior of the finite-size 
(moment equations and Markovian) systems is comparable to the behavior of WC 
solutions or to the behavior of the finite-size moment equations, which differs at any 
finite scale from WC's behavior as shown in section |3~2| Indeed, we observed in figure 
Fig.[6]that the Hopf bifurcation of the deterministic WC system, which is a degenerate 
point of the moment equations, expands in the codimension two diagram into two 
branches of generic Hopf bifurcations, that imply the presence of cycles in the finite- 
size system. These Hopf bifurcations either unfold a singularity corresponding to the 
actual eigenvalues ±A of WC Hopf bifurcation or to twice this value by application 
of the lemma [T] and can be either stable or unstable. Therefore the boundaries of the 
region where the system is in a periodic regime are modified as soon as the system 
size is finite. The same phenomenon appears close to LPa, in which case the diagram 
unfolds into a Hopf bifurcation, a saddle node bifurcation, a fold of limit cycles and a 
saddle homoclinic bifurcation. We will observe that in both cases, the behavior of the 
Markovian system actually reflects the properties of the finite-size moment equations, 
and that quasicycles will appear in the neighborhood of WC's Hopf bifurcation. 

Let us first have a closer look to the behavior of the moment equations in the 
neighborhood of I\ — I\ h m the finite-size system. WC's Hopf bifurcation gives rise 
to two Hopf bifurcations whose codimension two curves are tangent to the axis n = 0. 
The curve HI is related to a singular eigenvalue converging as N goes to infinity to the 
value of WC Hopf bifurcation, and gives rise to stable limit cycles with period close to 
the expected cycle in WC system. The other branch, H3, corresponds to faster cycles 
with period around half WC cycles' period. Let us now fix a certain size of the network 
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TV. In the moment equations, and let I\ be a free bifurcation parameter. As I\ increases, 
cycles will appear corresponding to the Hopf bifurcation HI at a value Ihi(N). We 
observe that this value is increasing as the network size increases. Therefore, in regions 
where WC system has a stable focus, the finite-size moment equations show oscillations 
of frequency close to WC cycles' frequency. For instance for n = 0.02 corresponding to 
the bifurcation diagram of figure Fig. [5] we have 7#i(50) = —3.37, defining a region 
of quasicycles I\ 6 [—3.37,-3.245]. This size of the zone where quasicycle appear 
decreases very fast (as we can see from the shape of the graph of HI in the codimension 
two plot Fig. [6] and for TV = 1000 we have Ij/i(50) = —3.248. It is important to note 
that at a size TV = 50, oscillations are very difficult to see because of the small system 
size. Moreover, in the particular case treated, the graph of the HI bifurcation is almost 
a vertical line: there is almost no difference between the value of HI across different 
system sizes. The same phenomenon is observed in the simulations: quasicycles close 
to the Hopf bifurcation point barely exist, and the value of I\ needs to be very close 
from 1 1 h to start observing these quasicycles. 

WC's family of limit cycle lose stability through a saddle-node homoclinic bifurca- 
tion. We observe that unfolding this bifurcation, we find out that the branch of limit 
cycles in the finite-size system undergo a fold of limit cycles corresponding to a decreas- 
ing value of I\ , having the effect of reducing the amplitude of the interval of I\ values 
related to a periodic orbit. The curve of fold of limit cycles is no longer vertical in that 
case and a wide non-trivial region distinguishes the finite-size and WC's periodic orbit 
parameter regions. We observe again in that case that the LPC2 line actually corre- 
sponds to a good estimate of the disappearance of cycles. Indeed, the codimension 2 
bifurcation diagram shows that for I\ = oscillations should disappear at a system 
size close to 550 neurons. In figure Fig. [IT] we observe that indeed, for TV" = 500 neurons 
and I\ = 0, the system does not show any oscillations, and for TV = 700 neurons the 
system starts showing an oscillatory behavior. Moreover, in that case we observe that 
the amplitude of the cycles is not small, consistent with the cycles we observe in the 
finite-size system. 

The presence of a chaotic orbit as seen in Fig. [7] and the fact that indeed in such 
regions the behavior of the system is aperiodic suggests that in these regions close to a 
limit cycle of Wilson and Cowan system, the underlying Markovian system can present 
a low-dimensional chaotic macroscopic behavior arising through the interaction of a 
large but finite number of neurons, a phenomenon that could be related to the work 
of El Boustani and Destexhe in [23] , where intrinsically stochastic neurons are shown 
to yield low-dimensional chaotic macroscopic behaviors. Results are presented in figure 
Fig. (IT] 

This particularly good correspondence between cycles of the moment equations 
and quasiperiodic orbits of the Markov chain is a very promising result. However, it 
is important to note that the additional cycles corresponding to the branches H2 and 
H4 do not appear in the Markov simulations. This might be linked with the stability 
of such orbits and the system size at which stable periodic orbit exist. This study 
therefore indicates that further exploration is needed to assess more clearly a definitive 
relationship between quasicycles and actual orbits of the moment equations, and there 
is a good chance that these two phenomena are related. This precise study is not within 
the scope of the present paper and we are currently actively working on this topic. Of 
high interest also is the fact that in some regions corresponding to the vicinity of a 
Hopf bifurcation, the moment equations present a low-dimensional chaotic behavior. 
Such phenomena, of interest in neuroscience, have, to the best of our knowledge, not 
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Fig. 11 Simulations of the Markovian system and periodic obits. Simulations are performed 
using Doob-Gillespie algorithm in order to simulate a large number of sample paths of the 
Markovian system, and the power spectrum presented are averaged across the realizations. 
For I2 = —5, orbits of the WC system exist for I± £ [—3.245,0.54], and around this zone, 
depending on the system size, cycles might exist or not (see text). 



been shown, and we believe that if these phenomena indeed appear in the Markovian 
system, this is of crucial interest in neuroscience and in the understanding of large-size 
stochastic systems. 
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4 Discussion 

In this paper, we studied the different solutions of newly developed models for describ- 
ing the mean-field limit of Markovian neural networks modeling large-scale cortical ar- 
eas: BCC and Bressloff models [TTlllO] . These equations were originally derived from a 
Markovian implementation of the firing of each neurons in the population, whose transi- 
tion probability satisfied a differential equation, the master equation of the process[22]. 
Both models studied drew on this framework but, since they addressed different neu- 
ronal regimes, dealt with different variables and different scaling, and after truncation 
of the moment hierarchy yielded sets two coupled equations, one governing the mean 
firing activity and the other governing the same-time correlations of the firing activity 
in the network (or the ordered cumulant that can be deduced from the correlation and 
the mean firing rate). The particular aim of both models was to recover, in the limit 
where the number of neurons is infinite, the standard WC equations. 

We transformed these models in order to study the proportion of active neurons 
in each populations, in networks composed of different populations of neurons with 
different dynamics and composed of different number of neurons. Both models include 
finite-size corrections, and the size of the network appears simply in the equations as a 
parameter of the model. The limit when the number of neurons tend to infinity appears 
in that view as a regular perturbation of the finite-size system in our particular choice 
of variables, and that both the BCC- and Bressloff-derived models converged towards 
the same equations as the number of neurons tends to infinity, the infinite-size system. 

We analyzed the solutions of the infinite-size system, and established the fact that 
all solutions of WC system were associated with solutions of the infinite-size system 
with null correlations. However, we showed that the stability of the solutions obtained 
in that fashion did not necessarily have the same stability properties as in WC system: 
fixed points of the Wilson and Cowan system correspond to fixed points of the infinite- 
size system with zero correlations having the same stability, but all the cycles of the 
WC system are either unstable or neutrally stable in the infinite-size system. This 
destabilization of cycles gives rise to new behaviors (fixed points or periodic orbits) 
with non-zero correlations, the correlation-induced behaviors. 

When studying the finite-size equations and the way it relates to the properties of 
the infinite-size system, we presented evidence that the mean-field limit appears to be 
a singular limit, in the sense that for any finite-size of the network, the system presents 
behaviors the infinite-size system does not feature. Moreover, singular points appear 
in the limit of infinitely many neurons, and for infinitely many neurons, branches of 
bifurcations collapse and disappear. Besides these qualitative distinctions between fi- 
nite and infinite-size systems, we numerically found mesoscopic-scale effects, namely 
behaviors that are only present at a certain population scale, and disappear for smaller 
or larger populations. The scale at which these behaviors took place coincided with the 
scale of typical cortical columns. Such behaviors are neither captured by the study of 
small-networks nor by the mean-field limit, stressing the importance of precise descrip- 
tions of these intermediate scales in tractable approaches. We want to emphasize that 
the system-size expansion and the truncations are not valid at the bifurcation points. 
However, the bifurcation analysis we performed allowed us to uncover the number and 
stability of fixed points and limit cycles corresponding to different parameters. The 
actual transition of the stochastic process is a very complex phenomenon which would 
be worth studying but that present very important technical intricacies. We also note 
that the underlying Markovian models are stochastic models. If the moment equations 
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present multiple attractors, the stochastic system will randomly switch between the 
different attraction basins. A large-deviation study in the case where there is a small 
parameter (Bressloff case and BCC case in the fully connected system treated in the 
manuscript) can allow one to identify the attractor switching and the mean exit time 
of each attraction basin in order to more precisely study the behavior of the underlying 
stochastic system, that could possibly be performed using Freidlin and Wentzell the- 
ory in the large N case |28| . As already mentioned, Paul Bressloff opened the way to 
such studies by using WKB approximations allowing to compute the transition rates 
between metastable attractors in a one-population model [llj . Alternative approaches 
could make use of the stochastic dynamical systems theory fBTJ33| . These approaches 
would allow directly studying the bifurcations of the initial stochastic system in the 
Langevin approximation. 

Correlation-induced and finite-size behaviors were then compared to the initial 
Markovian model. We observed that when WC system presented an attractive fixed- 
point behavior, each sample path of the Markov chain presented stochastic variations 
around these fixed points, the correlations tended towards zero and the mean coincided 
with WC fixed point. In that case, the network regime is therefore an asynchronous 
state, and is well approximated by Wilson and Cowan system. The validity of the 
moment equations close to bifurcations was shown in the derivations to be mainly valid 
far away from bifurcations. However, we showed that the codimension two bifurcations 
of the moment equations was closely matched with similar phenomena (e.g. bistability) 
observed in simulations of the Markovian system. However, we also observed that both 
BCC and Bressloff models can present additional behaviors in these regimes that are 
not present in the simulations of the Markov chain. The moment expansion breaks 
down when considering periodic orbits of the WC system, corresponding to quasi- 
periodic stochastic solutions of the stochastic system, and different techniques such as 
the stochastic phase reduction, Langevin approximation and Floquet theory are very 
promising. However, we provided numerical evidences showing that actual cycles in the 
moment equations could be related to the phenomenon of quasicycles in the Markovian 
system. These numerical observations strongly suggest further exploration along these 
lines. 

The precise study of all these points will further allow understanding the behavior 
of large scale systems, and will particularly be interesting in the study of rhythms in the 
brain, and possible synchronization or desynchronization as a function of the network 
size and of the connectivity. A particularly interesting point of the models studied in 
this paper is that they describe a stochastic process by a set of differential equations, 
where the size of the network appears as a perturbation parameter of a set of ordinary 
differential equations. This approach can be promising in order to study bifurcations 
of a stochastic processes and qualitatively analyzing its behaviors and the dependency 
of its solutions as a function of the parameters. We also noted that the behavior of 
each sample path of the process was possibly very different from the behavior of its 
first moments, which raises the question of the best variable to describe such stochastic 
processes. 

The mathematical method developed to study the stability of the correlation equa- 
tions, involving Kronecker products and sum and the application of Floquet theory 
is of potential importance in the study of networks with correlations, since the form 
of the correlation equation is quite standard in this kind of problems (see e.g. [35] in 
the case of small noise expansion). This approach can also be applied to the study of 
the stability of solutions in networks involving Hebbian learning, since the connectivity 
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weights evolution depends on the firing correlations and the stability of the weights 
involves very similar equations as the ones we studied here. 

The need for extensions of the present study to more realistic Markovian models 
involving more relevant biologically realistic transitions is a direct consequence of the 
present work. Indeed, we have seen that the transitions chosen in the Markovian model 
are only one-step jumps. In biological terms, this assumption means that when rij(t) 
neurons are active in population at time i, at least rii(t) — 1 neurons are active at time 
t + dt, i.e. keep firing, which is not biologically realistic since after firing, a neuron 
immediately returns to quiescence. Therefore extending the framework to a more com- 
plex state space of the Markov chain would be of great interest for the development 
of these models towards more realistic biological grounds, and would probably allow a 
better assessment of the properties of cortical areas. 

Finally, in order to go beyond the description of the moments of the Markov chain, 
we believe that one shall turn towards alternative ways to derive mean-field equations 
for such systems. A particularly interesting way to study such systems we are currently 
exploring builds upon recent methods developed for the study of Markov chains, mainly 
in the field of queueing theory, such as fluid limits and large deviations approaches. The 
fluid limit technique is based on a space-time proper rescaling depending on the process 
under consideration, that yields deterministic or stochastic limits of the activity of the 
network when the number of neurons tends to infinity. Large deviation techniques allow 
one to derive the finite-size corrections (see e.g. [52] in a case of stochastic ion channels) 
and is a natural extension of our work. 

Acknowledgments: This work was supported by NSF DMS 0817131. 



A An alternative derivation of the moment equations 

In this appendix we show that the moment equation corresponding to the Markov chain gov- 
erned by equations |2jl can be derived from a Rodriguez-Tuckwell moment expansion on the 
Langevin approximation of the Markov chain. Following Kurtz approach in 1341 . we know that 
the dynamics of the our rescaled variables pi = m/Ni, which is a Markov chain in the initial 
setting, approaches as N goes to infinity a continuous diffusion process (Xt)t>o satisfying the 
equation: 

dX\ = ( - at X\ + fi(si(t,X t ))) dt + ^ a, XI + nis.it, X t ))dW l t (13) 

where S,(t,X t ) = £jLi w^X* + hit). 

Following the works of Rodriguez and Tuckwell [44145] . we derive from this equation the 
dynamical system governing the approximate moments of X. Denoting by rrij the mean value 
of Xj and by dj the correlation between Xi and Xj , a direct application of Rodriguez and 
Tuckwell formula applied to our particular form of dynamical system yields: 

-^i = -ajnij + fj ( Sj ) + |/j'(sj) J2f=i J2p=i w ]pWjl C lv 
, *°SL = -( ai + a] )C l3 + f({ Si ) E£i WuCy + f'fa) ElLi WjiCu 

Where s; = T*lj— i w ij m t + h{t). These equations therefore appear as a perturbation of the 
BCC equations with an additional nonlinear term in the dynamics or order two. Truncation at 
order 1 in the small parameter l/N yields exactly BCC equations. Moreover, it is interesting 
to note that these equations again correspond in the limit N — > oo, to the infinite model 
described in section [2] 
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B Kronecker Algebra: Some useful properties 

In this appendix, we review and prove some useful properties of Kronecker products of matrixes. 
We recall the definition of the the function Vect transforming a M X JV matrix into a M N- 
dimensional column vector, as defined in 141 1 : 



Vect 



R MxN ^ R MN 

X k> [Xj X M1 , Xia (t) , . . . , X M2 (t) , ... X 1Ar (t) X MN (t)] 3 



Let us now denote by © the Kronecker product defined for A £ R mxn and B £ R rxs as the 
(mr) X (us) matrix: 

/ onB ai 2 S ■ ■ ■ a ln B \ 
a2iB a22^ ■■■ a,2 n B 



A®B 



V o-miB a m 2B ■ ■ ■ a mn B ) 



For standard definitions and identities in the field of Kronecker products, the reader is referred 
to [12]. We recalled in the main text the following identities for A,B,D,G,X e R MxM , I M 
be the M X M identity matrix and A ■ B or AB denote the standard matrix product: 

(Vect(AXB) = (B T <g> A)Vect(X) 
A®A = A®I m +Im®A . (14) 

(A ® B) ■ (D ® G) = (A ■ D) x (B ■ G) 

The relationship © is called Kronecker sum. 

Proposition 2 Let A and B in R MxM , and assume that A has the eigenvalues {A;; i = 
1, . . . , M } and B the eigenvalues i = 1, , M}. Then we have: 

— A® B has the eigenvalues (i, j) £ {1, . . . , M} 2 }. 

— A (B B has the eigenvalues {\i + fij; (i, j) 6 {1, . . . , M} 2 }. 

Proof Let Ai (resp./ij) be an eigenvalue of A (resp. B) with eigenvector u (resp. v), and define 
the matrix z = u ■ v T (i.e. mu = «i iij). We have: 

© Vect(z) = Vect (A ■ u ■ v T + u ■ v T ■ B T ^j 

= Vect(A,u -v T +u-(B ■ v) T ^ 

= AiVect(z) + Vect(u(n j v T ) S ) 
= (Ai + /ij)Vect(z). 

which entails that Vect(z) is an eigenvector of A(u) © A(v) associated with the eigenvalue 
Ai + \j . Similarly, we have for z = v ■ u T : 

(A © Vect(z) = Vect^B ■ v ■ u T A T ) 

= Vect((Ayv) ■ (Ai u T )) 
= Ai^ h Vect(z) 

The dimension of A © B and A © B is M 2 and there are exactly M 2 linearly independent 
matrices z possible built in the proposed fashion, therefore we identified all possible eigenvalues. 

Theorem 3 Let 'Pit) be the solution of the matrix differential equation: 

\#(0) = I M 
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and the solution of: 

=(A(t)®A(t)Wt) 



*(0) = I M 2 



We have <P(t) = <P(t) ® &(t). 



Proof Indeed, we have, using the basic properties of the Kronecker product recalled in equation 
§ 

d<P(t) ® $(t) d<P(t) , , d<P(t) 

= (A{u(t)) ■ 0(1)) Cg 0(t) + 0(f) Cg) (A(i/(t)) ■ 0(f)) 

= {A(v(t)) ■ 0(f)) Cg (7 M ■ *(*)) + (/M ■ *(*)) ® (A(l/(t)) ■ 0(f)) 

= (A(i/(t)) Cg 7 M ) ■ (0(f) ® 0(t)) + (7m ® A(i/(t))) ■ (0(f) Cg 0(f)) 

= (A(i/(t)) eg 7 M + l M eg A(i/(t))) ■ (0(f) eg 0(f)) 

= (i(Kt))9 4Mt)))-(#(t)®*(i)) 

Therefore 0(t)cg0(t) satisfies the same differential equation as •P'(f) and moreover, 0(O)(g0(O) = 
Im®Im = 7 M 2, and therefore by existence and uniqueness of the resolvent, •f'(t) = 0(f)® 0(f). 



C Genericity of the Hopf bifurcation found 

In this appendix we derive the expression of the first Lyapunov exponent of the bifurcation, 
which proves that the existence of the Hopf bifurcation exhibited in section |3.1.2| In that 
section, we derived the expression of the Jacobian matrix at the considered fixed point: 

-a + wf'(s) \f"{s)w 2 
2f'(s)% 2(-a + w/'W) 

At the bifurcation point, we have —a* +«)/'(«) = 0, and therefore at this point the Jacobian 
matrix reads: 

\f"{I)w^ 
.2/'(7)f 



To 



The eigenvalues of this matrix under the assumptions of section |3.1.2| arc iiairj where u>o 
y/—f'(I)f"(I)w 3 /N. We define q the right eigenvector of Jo associated with icoo: 



s/2f'{I)w/N 

and p the right eigenvector of Jq associated with the eigenvalue iu>o: 

l) 



For the sake of simplicity, we also name the components of the vector field of the system: 

/!((£),«*) =wf'{I)nu + f{wv + I) + \f"(wv + I)w 2 C 
/ 2 ((£),a*) = -2wf'(I)C + 2f'(wu + I)w(C+%) 

Following L 35], we define #(Q), and C(Q), (**), (*|)) the second and third deriva- 

tives of the vector field, which are bi- and tri-lincar forms. We have the following expressions 
for these multilinear functions: 
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{ Bi O • O) = w2 /" w xi 3:2 + 1 /"W ™ 2 («* w + » 1 X2 ) 

I B2(C;)-©) = 2/"(/) «> 2 (xj J/2 + 2/1 xa) + £ /"(/) ™ 2 xi x 2 

I Cl ( («) > (w) ' (w) ) = / <3) W W3 *l*a»3 + | /< 4 > (/) 1« 4 (XI Z 2 J/3 + XI y 2 X3 + 2/1 X 2 X 3 ) 
[C2((^), (£), (II)) = jr f (3) (I) «> 3 X! x 2 x 3 + 2 /< 3 )(/) ™ 3 (xi x 2 y 3 + X! 2/ 2 x 3 + W x 2 x 3 ) 

We are now in position to compute the first Lyapunov exponent (i(0) using the formula: 
h(0) = Rc((p, C(<?, q,q)) - 2(p, B(q, J^ 1 B(q,q))) + (p, B(q, (2 i £j /d - J ) _1 B( 9) </))>) 

2w 

where (x, J/} denotes the complex inner product x T • 1/ and the sum of three terms denoted A, 
B and C are the real parts of the terms involved in the expression of the Lyapunov exponent. 
After straightforward but tedious calculations (that can be conveniently performed using a 
formal calculation tool such as Maple), we obtain: 

« _ w 2 N ( f"(lf l / (3) (//'(-T)) , ?fH(r\A 

a - 7W V w 2 wo +2/ V) ) 

r ™ 2 N . 2 /"( J ) 2 

TUT 1 3 

which yields the expression for the Lyapunov exponent: 

Zj(0) = — [A-2B + C] 
2uiq 

f (s) (l)f'{i) +s"( I ? (— - ^f) 

V wo/ V^o 3 /. 



D Finite-size effects in BCC two-populations Model I 

In this appendix, we study the two-populations BCC system corresponding to Model I and 
show that the finite-size effects are closely related to what is observed in Bressloff model as 
studied in section [3.2| BCC finite-size equations read: 

' a\ = -aai + /(si) + § f"(si) (w 2 L1 en + wf 2 c 2 2 + 2w 12 wn c 12 ) 
a' 2 = -aa 2 + f(s 2 ) + | /"(s 2 ) {w\ 2 c 22 + w\ x en + 2 11)22 IU21 C12) 
( 4i = -2ocn + 2/'(si) (ton en + tt)i2 C12) + 2/'(si) ai ra e 
c 22 = -20C22 + 2f'(s 2 ) (W21 C12 + «>22 C22) + 2 n; /'(S2) W22 a 2 
c' 12 = -2oci2 + /'(si) (Wll C12 + 1U12 C22 + a2 W12 «i) ■ • ■ 
+/'(«s) (">21 en + 1022 C12 + W21 ai n e ) 

where we denoted: 

si = u>n ai + UI12 a 2 + ii 
S2 = W21 a, 2 + W22 a 2 + 12 

Similarly to Bressloff case, BCC model features two families of limit cycles. One of these 
branches corresponds exactly to the branch of limit cycles of WC system starting from a 
Hopf bifurcation and disappearing through a homoclinic bifurcations. Two additional Hopf 
bifurcations appear related to the family of periodic orbit corresponding to the correlation- 
induced cycle evidence in the analysis of the infinite-size system. This branch of limit cycles 
exist whatever n, and loses stability through a Neimark-Sacker bifurcation as the number of 
neurons increases. This bifurcation generates chaos for large enough networks, which clearly 
does not exists in WC system. 
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Fig. 12 Bifurcation diagram for the BCC system. Blue lines represent the equilibria, pink 
lines the extremal values of the cycles in the system. Bifurcations of equilibria are denoted 
with a red star, LP represents a saddle- node bifurcation (Limit Point), H a Hopf bifurcation. 
The four Hopf bifurcations share two families of limit cycles. The branch corresponding to the 
smaller values of ii undergoes two fold of limit cycles, and the other branch of limit cycle a 
Neimark Sacker (Torus) bifurcation. 



This family of limit cycles presents stability for networks containing between 100 and 
17.000 neurons, corresponding to a finite size effect that appears only in the region of interest 
of the cortical columns. As the number of neurons increase, the system keep the same number 
and type of bifurcations, and the infinite model appears to be a singular case where different 
bifurcations meet (see Fig. |13| , and very large networks lose the property of presenting a stable 
cycle, and present a chaotic behavior. 
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